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September  14,  201 2 


Re:  Contract  No.  W91 1NF-12-I-0144  -  award  to  Colorado  State  University  (proposal  61946-EL-CF)  to 
fund  $7,000.00  for  Conference  Support  Grant  for  FEM2012  International  Workshop  (Dr.  Dev  Palmer. 
Program  Manager) 


Dear  Colleagues: 

Please  find  enclosed  three  copes  of  the  Conference  Proceedings  of  FEM2012. 

Please  also  note  that  the  signed  Final  Certification  statement  is  being  sent  to  ARC)  as  well,  at  the  address 
indicated  in  the  contract. 


The  conference  support  grant  in  the  amount  of  $7000  by  ARO  was  used  to  support  the  conference  student 
paper  program  that  was  established  for  this  conference  for  the  first  time.  The  ten  students  winners  of  the 
FEM2012  Student  Paper  Competition  received  $700  student  travel  support  grants.  The  student  awardees 
were  recognized  at  the  conference  banquet  (each  of  the  students  received  a  plaque  at  the  banquet),  and 
presented  their  winning  papers  at  the  workshop. 

I  hope  that  upon  the  receipt  of  the  copies  of  the  Conference  Proceedings  and  the  signed  Final  Certification 
statement  the  ARO  will  be  able  to  release  the  final  payment  of  $700  of  the  grant  (as  stipulated  in  the 
contract). 

Thank  you  very  much  for  the  ARO's  support  of  our  conference  and  student  program. 

Sincerely, 


General  Chair  of  the  1 1th  International  Workshop  on  Finite  Elements  for  Microwave  Engineering - 
FBM2012,  held  from  June  4-6,  2012,  in  Estes  Park,  Colorado 
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Dr.  Branislav  M.  Notaros.  Professor 

Department  of  Electrical  and  Computer  Engineering 

Colorado  State  University 

1373  Campus  Delivery 

Fort  Collins,  CO  80523-1373 

Phone;  (970)  491-3537 

Fax:  (970)  491-2249 

E-mail:  notarns@.colostate.ed u 

Web:  w  n  u  cnerxolostate.edu  -notaros 
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U,S.  Army  Contracting  Command  -  Aberdeen  Proving  Ground 
Research  Triangle  Park  Division 
ATTN:  CCAP-SCR  (CLOSEOUTS) 

P.O.  Box  12211 

Research  Triangle  Park,  NC  27709-221 1 

Re:  Contract  No,  W91  INF-12-1-0144  -  award  to  Colorado  State  University  (proposal  61946-EL-CF)  to 
fund  $7,000,00  for  Conference  Support  Grant  for  FEM2012  International  Workshop  (Dr,  Dev  Palmer. 
Program  Manager) 

Dear  Colleagues: 

This  is  the  signed  Final  Certification  statement  for  the  above  ARO  Grant  for  FFM201 2  Conference. 

Please  also  note  that  the  three  copes  of  the  Conference  Proceedings  of  FEM2012  were  sent  to  ARO  as 
well,  at  the  address  indicated  in  the  contract. 

The  conference  support  grant  in  the  amount  of  $7000  by  ARO  was  used  to  support  the  conference  student 
paper  program  that  was  established  for  this  conference  for  the  first  time.  The  ten  students  -  winners  of  the 
FEM2012  Student  Paper  Competition  received  $700  student  travel  support  grants.  The  student  awardees 
were  recognized  at  the  conference  banquet  (each  of  the  students  received  a  plaque  at  the  banquet),  and 
presented  their  winning  papers  at  the  workshop. 

1  hope  that  upon  the  receipt  of  the  copies  of  the  Conference  Proceedings  and  the  signed  Final  Certification 
statement,  the  ARO  will  be  able  to  release  the  final  payment  of  $700  of  the  grant  (as  stipulated  in  the 
contract). 

Thank  you  very  much  for  the  ARO*  s  support  of  our  conference  and  student  program. 


Final  Certification  statement: 

/  certify  that  the  funds  provided  by  the  Army  Research  Office  under  Grant  W9 11  NF-i  2-1-01 44  were 
expended  in  accordance  with  the  provisions  of  the  grant  and  that  required  conference  proceedings 
have  been  delivered  to  the  Army  Research  Office. 


Branislav  Notaros 


General  Chair  of  the  1 1th  International  Workshop  on  Finite  Elements  for  Microwave  Engineering  - 
FEM2012.  held  from  June  4-6.  20 12,  in  Estes  Park.  Colorado 
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Department  of  Electrical  and  Computer  Engineering 

Colorado  State  University 

1373  Campus  Delivery 

Fort  Collins,  CO  80523-1373 
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Fax:  (970)  491-2249 
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Web:  www.engr.colostate.edu/~notaros 


For  three  copes  of  the  Conference  Proceedings  of  FEM2012 


United  States  Postal  Service 
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P.O.  Box  12211 
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Physical  Address 


U.S.  Army  Research  Office 
ATTN:  RDRL-ROE-L 
4300  South  Miami  Boulevard 
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Main  Phone  Number:  (919)  549-0641 
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U.S.  Army  Contracting  Command  -  Aberdeen  Proving  Ground 
Research  Triangle  Park  Division 
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FEM2012 


Student  Paper  Award  and  US  Army 
Research  Office  Travel  Support  Grant 

is  presented  to 

Thomas  Bauernfeind 

for  a  paper  at 

The  11th  International  Workshop  on  Finite 
Elements  for  Microwave  Engineering 

June  4-6,  2012,  Estes  Park,  Colorado,  USA 


11th  International  Workshop  on  Finite  Elements 
for  Microwave  Engineering  -  FEM2012 

June  4-6,  2012,  Estes  Park,  Colorado,  USA 


http://www.engr.colostate.edu/FEM2012 
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Technical  and  Social  Program  at  a  Glance 


Sunday,  June  3,  2012 

Welcome  Reception  19:00-21:00  (Music  Room) 


Monday,  June  4,  2012 

Breakfast  7:00-8:00  (MacGregor  Room) 

Technical  Sessions  8:00-11:50  (coffee  break  9:40-10:10) 

Lunch  11:50-13:10  (MacGregor  Room) 

Technical  Sessions  13:10-17:20  (coffee  break  15:10-15:40) 
Stanley  Hotel  Ghost  &  History  Tour  -  Monday  Tour  17:30-19:00 
Conference  Banquet  19:30-22:30  (MacGregor  Room) 


Tuesday,  June  5,  2012 

Breakfast  7:00-8:00  (MacGregor  Room) 

Technical  Sessions  8:00-11:50  (coffee  break  9:40-10:10) 

Lunch  11:50-13:10  (MacGregor  Room) 

Technical  Sessions  13:10-17:00  (coffee  break  15:10-15:40) 
Stanley  Hotel  Ghost  &  History  Tour  -  Tuesday  Tour  17:30-19:00 


Wednesday,  June  6,  2012 

Breakfast  7:00-8:00  (MacGregor  Room) 

Trip  to  Rocky  Mountain  National  Park  8:00-14:00  (picnic  lunch) 
Open  Forum  Discussions  15:00-18:00  (coffee  break  16:30-16:50) 
Closing  Reception  18:00-21:00  (The  Elkhorn  Lodge  &  Ranch) 


The  FEM  International  Workshop  is  a  highly  focused  biannual  event  providing  an  ideal 
meeting  place  for  researchers  who  are  active  in  either  the  theoretical  and  numerical  development 
of  finite  element  and  hybrid  methods  or  their  application  to  a  broad  range  of  electromagnetic 
problems. 

The  FEM2012  Technical  Program  combines  82  excellent  papers  organized  in  12  special 
sessions.  It  also  features  three  open  forum  discussions,  moderated  by  leading  FEM  researchers 
and  practitioners  from  academia,  industry,  and  government. 

The  FEM2012  Social  Program  includes  welcome  reception,  conference  banquet,  closing 
reception,  trip  to  Rocky  Mountain  National  Park,  and  Stanley  Hotel  Ghost  &  History  Tour. 
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FEM2012  Committees 

General  Chairs: 

Branislav  M.  Notaros,  Colorado  State  University 
notaros@colostate.edu 

Jian-Ming  Jin,  University  of  Illinois 
i-iinl@ad.uiuc.edu 


Technical  Program  Committee  Chair: 

Milan  M.  Ilic,  University  of  Belgrade  and  Colorado  State  University 
milanilic@,etf.rs 


Local  Organizing  Committee  Chair: 

Olivera  Notaros,  Colorado  State  University 
olivera@colostate.edu 


Scientific  Committee: 

A.  C.  Cangellaris,  University  of  Illinois  Urbana-Champaign,  USA 
Z.  J.  Cendes,  Ansoft  LLC.,  USA 

D.  B.  Davidson,  University  of  Stellenbosch,  South  Africa 
R.  Dyczij-Edlinger,  Universitat  des  Saarlandes,  Germany 
R.  L.  Ferrari,  Trinity  College,  UK 
J.  Jin,  University  of  Illinois,  Urbana-Champaign,  USA 

L.  Kempel,  Michigan  State  University,  USA 
J.-F.  Lee,  Ohio  State  University,  USA 

R.  Lee,  Ohio  State  University,  USA 

G.  Pelosi,  University  of  Florence,  Italy 

T.  Rylander,  Chalmers  University  of  Technology,  Sweden 

M.  Salazar-Palma,  Universidad  Carlos  III  de  Madrid,  Spain 
J.  L.  Volakis,  Ohio  State  University,  USA 

M.  N.  Vouvakis,  University  of  Massachusetts,  USA 

J.  P.  Webb,  McGill  University,  Canada 

T.  Weiland,  Technische  Universitaet  Darmstadt,  Germany 


2 


FEM2012  Financial  Sponsors 


ANSYS  (ANSOFT) 


ANSYS 


CST  -  Computer  Simulation  Technology 


EM  Software  &  Systems-S.A.  (Pty)  Ltd.,  FEKO 


Comprehensive  Electromagnetic  Solutions 


FEM2012  Student  Paper  Competition 
and  Travel  Grants  Sponsor 


US  Army  Research  Office  (ARO) 


The  views,  opinions,  and/or  findings  contained  in  this  publication  are  those  of  the  authors  and  should 
not  be  construed  as  an  official  Department  of  the  Army  position,  policy,  or  decision,  unless  so 
designated  by  other  documentation. 


FEM2012  Technical  Sponsors 

IEEE  (Institute  of  Electrical  and  Electronics  Engineers) 
IEEE  Antennas  and  Propagation  Society 
IEEE  Microwave  Theory  and  Techniques  Society 
University  of  Illinois 

URSI  (International  Union  of  Radio  Science) 


m 

MTT-S" 
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FEM2012  Student  Paper  Competition  Awardees 


These  ten  students  are  winners  of  the  FEM2012  Student  Paper  Competition  and  will  receive 
$700  student  travel  support  grants  sponsored  by  the  US  Army  Research  Office  (ARO).  The 
student  awardees  will  be  recognized  at  the  conference  banquet,  and  must  present  their  winning 
paper  at  the  workshop. 


Student  name 

Affiliation 

Paper  number  and  title 

Thomas  Bauemfeind 

Institute  for  Fundamentals  and  Theory 
in  Electrical  Engineering,  Graz 
University  of  Technology,  Inffeldgasse 
18,  8010  Graz,  Austria 

1.6.  Investigation  of  the  Scattering 

Behavior  of  Transponder  Antennas  in 

Case  of  Nonlinear  Termination 

Stylianos  Dosopoulos 

Electroscience  Laboratory,  Ohio  State 
University,  1320  Kinnear  rd., 

Columbus  OH  43212,  USA 

2.3.  Non-Conformal  and  Parallel  Interior 
Penalty  Discontinuous  Galerkin  Method 
for  the  Time-Domain  Maxwell’s 

Equations 

Grzegorz  Fotyga 

Gdansk  University  of  Technology 
Department  of  Microwave  and  Antenna 
Engineering  11/12  Narutowicza  Street, 
80-233  Gdansk,  Poland 

5.2.  Efficient  Model  Order  Reduction  for 
FEM  Analysis  of  Waveguide  Structures 
and  Resonators 

Qing  He 

School  of  Electrical  and  Computer 
Engineering,  Purdue  University,  465 
Northwestern  Avenue,  West  Lafayette, 
IN  47907,  USA 

3.5.  An  Explicit  and  Unconditionally 

Stable  Time-Domain  Finite-Element 
Method  of  Linear  Complexity 

Gergely  Koczka 

Institute  for  Fundamentals  and  Theory 
in  Electrical  Engineering,  Graz 
University  of  Technology,  Inffeldgasse 
18,  8010  Graz,  Austria 

6.8.  An  Iterative  Domain  Decomposition 
Method  for  Solving  Wave  Propagation 
Problems 

Ana  Manic 

Colorado  State  University,  Electrical  & 
Computer  Engineering  Department, 

1373  Campus  Delivery,  Fort  Collins, 

CO  80523-1373,  USA 

9.7.  Symmetric  and  Nonsymmetric  FEM- 
MoM  Techniques  Using  Higher  Order 
Hierarchical  Vector  Basis  Functions  and 
Curved  Parametric  Elements 

Safa  Salman 

Electroscience  Laboratory,  Ohio  State 
University,  1320  Kinnear  rd., 

Columbus  OH  43212,  USA 

1 . 1  .Highly  Coupled  Bowtie  Array  Design 
with  Passive  Elements 

Nada  Sekeljic 

Colorado  State  University,  Electrical  & 
Computer  Engineering  Department, 

1373  Campus  Delivery,  Fort  Collins, 

CO  80523-1373,  USA 

8.3.  Rules  for  Adoption  of  Expansion  and 
Integration  Orders  in  FEM  Analysis 

Using  Higher  Order  Hierarchical  Bases 
on  Generalized  Hexahedral  Elements 

Ming-Feng  Xue 

Electromagnetics  Laboratory, 
Department  of  Electrical  and  Computer 
Engineering,  University  of  Illinois  at 
Urbana-Champaign,  Urbana,  IL  61801, 
USA 

6.7.  Preliminary  Research  on  the 
Discontinuous  Enrichment  Method  Based 
Domain  Decomposition  Scheme  for 

Solving  the  Three-Dimensional  Vector 
Curl-Curl  Equation 

Wang  Yao 

University  of  Illinois  at  Urbana- 
Champaign,  Electrical  and  Computer 
Engineering,  1406  West  Green  Street, 
Urbana,  IL  61801,  USA 

11.5.  Application  of  Tree-Cotree  Splitting 
to  the  Dual-Primal  Finite  Element 

Tearing  and  Interconnecting  Method  for 
Solving  Low-Frequency  Breakdown 
Problems 

The  views,  opinions,  and/or  findings  contained  in  this  publication  are  those  of  the  authors  and  should  not  be  construed  as  an 
official  Department  of  the  Army  position,  policy,  or  decision,  unless  so  designated  by  other  documentation. 
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Overall  Technical  Program 


The  Historic 
Stanley  Hotel 

Monday 
June  4 

Monday 
June  4 

Tuesday 
June  5 

Tuesday 
June  5 

Wednesday 
June  6 

Time 

Slots 

Billiard 

Room 

Music 

Room 

Billiard 

Room 

Music 

Room 

Music 

Room 

8:00 

Slot  1 

Paper  1.1 

Paper  2.1 

Paper  7.1 

Paper  8.1 

8:20 

Slot  2 

Paper  1.2 

Paper  2.2 

Paper  7.2 

Paper  8.2 

8:40 

Slot  3 

Paper  1.3 

Paper  2.3 

Paper  7.3 

Paper  8.3 

9:00 

Slot  4 

Paper  1.4 

Paper  2.4 

Paper  7.4 

Paper  8.4 

9:20 

Slot  5 

Paper  1.5 

Paper  2.5 

Paper  7.5 

Paper  8.5 

9:40 

Coffee  break 

10:10 

Slot  6 

Paper  1.6 

Paper  2.6 

Paper  9.2 

Paper  8.6 

10:30 

Slot  7 

Paper  1.7 

Paper  4.1 

Paper  9.3 

Paper  8.8 

10:50 

Slot  8 

Paper  3.1 

Paper  4.2 

Paper  9.4 

Paper  10.1 

11:10 

Slot  9 

Paper  3.2 

Paper  4.3 

Paper  9.5 

Paper  10.2 

11:30 

Slot  10 

Paper  3.3 

Paper  4.4 

Paper  9.6 

Paper  10.3 

11:50 

Lunch  break 

13:10 

Slot  1 1 

Paper  3.4 

Paper  4.5 

Paper  9.7 

Paper  10.4 

13:30 

Slot  12 

Paper  3.5 

Paper  4.6 

Paper  11.1 

Paper  12.1 

15:00-15:45 

13:50 

Slot  13 

Paper  3.6 

Paper  6.1 

Paper  11.2 

Paper  12.2 

Open  forum 

14:10 

Slot  14 

Paper  5.1 

Paper  6.2 

Paper  11.3 

Paper  12.3 

discussion  1 

14:30 

Slot  15 

Paper  5.2 

Paper  6.3 

Paper  11.4 

Paper  12.4 

15:45-16:30 

14:50 

Slot  16 

Paper  5.3 

Paper  6.4 

Paper  11.5 

Paper  12.5 

Open  forum 

15:10 

Coffee  break 

discussion  2 

15:40 

Slot  17 

Paper  5.4 

Paper  6.5 

Paper  11.6 

Paper  12.6 

16:30-16:50 

16:00 

Slot  18 

Paper  5.5 

Paper  6.6 

Paper  11.7 

Paper  12.7 

Coffee  break 

16:20 

Slot  19 

Paper  5.6 

Paper  6.7 

Paper  11.8 

Paper  12.8 

16:50-18:00 

16:40 

Slot  20 

Paper  5.7 

Paper  6.8 

Paper  11.9 

Paper  12.9 

Open  forum 

17:00 

Slot  21 

Paper  5.8 

Paper  6.9 

Committee  Meeting 

discussion  3 

17:20 

Technical  program  ends  for  the  day 

51  Modeling  and  Design  of  Antennas  and  Arrays  Using  FEM  (R.  Kindt,  D.  Filipovic,  E.  Topsakal) 

52  Discontinuous  Galerkin  Methods  (Q.  H.  Liu,  W.  Cai) 

53  Time-Domain  FEM  and  Applications  (D.  White,  O.  Biro,  T.  Rylander) 

54  FEM  and  CEM  Applications  in  Optics  and  Nanophotonics  (E.  Simsek) 

55  Model  Order  Reduction  Techniques  (R.  Dyczij-Edlinger,  V.  de  la  Rubia) 

56  Advances  in  Vector  Bases  for  CEM  (R.  Graglia,  A.  Peterson,  D.  Wilton) 

57  FETD  Modeling  of  Complex  Media  and  Structures  (R.  Lee,  F.  Teixeira) 

58  Adaptive  FEM,  Higher  Order  Bases,  and  Advanced  FEM  Formulations  (J.  Webb,  M.  Gavrilovic,  M.  Salazar  Palma,  L.  Garcia-Castillo) 

59  Advances  in  Hybrid  Methods  and  Multiphysics  Problems  (B.  Shanker,  L.  Kempel) 

510  FEM  Modeling  and  Applications  of  Metamaterials  and  Periodic  Media  (K.  Sertel,  J.  Volakis) 

511  Domain-Decomposition  Methods  (J.-F.  Lee,  Z.  Peng) 

512  Advanced  FEM/MoM  Modeling,  Design,  and  Optimization  (C.J.  Reddy,  U.  Jakobus,  J.  Zapata,  J.  Gil,  D.  Jiao,  A.  Cangellaris) 


Open  forum  discussion  1:  Open  Forum  Discussion  on  FEM  and  CEM  Bases,  Elements,  and  Formulations 
(A.  Peterson,  D.  Wilton,  J.  Webb,  K.  Sertel) 

Open  forum  discussion  2:  Open  Forum  Discussion  on  DDM,  MOR,  and  Efficient  FEM  and  Hybrid  Solutions 
(J.-F.  Lee,  Zhen  Peng,  R.  Dyczij-Edlinger,  V.  de  la  Rubia) 

Open  forum  discussion  3:  Open  Forum  Discussion  on  FEM-Based  Modeling,  Design,  and  Applications 
(L.  Kempel,  B.  Shanker,  R.  Kindt,  C.J.  Reddy,  K.  Zhao,  T.  Euler) 
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Presentation  Information 

Language 

The  working  language  of  the  workshop  is  English,  which  includes  all  presentations,  printed 
materials,  and  discussions. 

Presentation  Length  and  Format 

All  presentations  are  oral  and  are  allotted  20-minute  time  slots,  which  includes  15  minutes  for 
the  presentation  and  5  minutes  for  questions  and  discussion. 

All  presenters  should  contact  the  session  chairs  no  later  than  1 5  minutes  before  the  start  of  the 
session  to  confirm  their  presence  and  upload  their  presentation  files  to  the  session  computer.  If 
the  session  starts  immediately  after  the  end  of  the  previous  session,  presenters  should  arrive  to 
the  session  room  and  contact  the  session  chairs  during  the  coffee  break  or  lunch  break  preceding 
the  start  of  the  session. 

Technical  Equipment 

Each  session  room  is  equipped  with  a  laptop  computer  connected  to  a  video  projector.  Video 
resolution  is  1280x800  pixels.  Session  room  computers  run  under  Windows  7  Pro  and  have 
Microsoft  Office  2010  installed. 

File  Formats 

The  only  file  formats  supported  by  the  workshop  computers  are  Microsoft  PowerPoint  (*.ppt, 
*.pptx)  and  Adobe  PDF  (*.pdf). 

File  Upload 

Speakers  are  asked  to  upload  their  presentation  files  to  the  session  room  computer  no  later  than 
15  minutes  before  the  session  start.  If  the  session  starts  immediately  after  the  end  of  the  previous 
session,  speakers  should  upload  their  files  during  any  breakfast,  coffee,  or  lunch  break  preceding 
the  session  (or  even  the  day/evening  before). 

Session  room  computers  are  equipped  with  USB  ports  and  DVD  drives.  Staff  personnel  will  be 
available  to  assist  with  the  upload  procedure. 

Session  Chairs 

Special  session  organizers  will  serve  as  session  chairs.  Session  chairs  must  arrive  to  the  session 
room  no  later  than  15  minutes  before  the  session  start.  If  the  session  starts  immediately  after  the 
end  of  the  previous  session,  session  chairs  should  come  to  the  room  and  check  if  everything  is 
ready  during  the  coffee  break  or  lunch  break  preceding  the  start  of  the  session. 
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Conference  Floor  Plan  and  Hotel  Directions 


Sessions  are  held  in  Billiard  Room  and  Music  Room,  registration  and  coffee  breaks  are  in  Pinion 
Room,  breakfasts,  lunches,  and  conference  banquet  are  in  MacGregor  Room. 
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DRIVING  FROM  ALL  DIRECTIONS 
Take  Interstate  25  to  Exit  257B  toward  the  city  of 
Loveland  and  Estes  Park  onto  state  HWY  34.  Continue 
through  Loveland  into  the  Big  Thompson  Canyon.  As 
you  enter  Estes  Park,  turn  right  onto  Steamers  Parkway 
to  The  Stanley  gatehouse  entrance. 


The  Historic  Stanley  Hotel  is  located  just  over 
an  hour  northwest  of  Denver,  Colorado  in  the 
town  of  Estes  Park,  gateway  to  Rocky  Mountain 
National  Park. 


For 


detailed  directions  to  and  from  Denver  International 
Airport,  (DIA)  please  visit  our  website  at 
www.stanleyhotel.com. 


HOTEL 


333  East  Wonderview  Avenue  •  Estes  Park,  Colorado  80517 


(rj/e*tcr  /909 
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About  the  Venue 

The  Historic  Stanley  Hotel 


The  Stanley  Hotel,  http://www.stanleyhotel.com,  is  a  historic  landmark  hotel  in  a  spectacular 
mountain-view  location,  offering  old-world  charm  matched  with  the  latest  of  modern  amenities. 
Located  at  an  elevation  of  7,500  feet  (2,286  meters)  and  within  sight  of  the  Rocky  Mountain 
National  Park,  it  was  built  in  1907-1909  by  the  inventor  of  the  Stanley  Steamer  automobile 
(which  is  on  display  in  the  hotel  lobby).  It  inspired  Stephen  King  to  write  The  Shining,  and 
served  as  a  location  site  for  several  movies  (including  Dumb  and  Dumber).  It  hosted  many  US 
presidents,  royalties,  and  celebrities.  Many  believe  that  the  hotel  is  haunted  (the  most  cases  of 
ghostly  activity  were  reported  for  the  ballroom  in  which  we  will  have  our  banquet  dinner).  The 
popular  Stanley  Hotel  Ghost  &  History  Tour  is  included  in  the  workshop  package. 


Estes  Park 

Estes  Park,  http://www.estesparkcvb.com,  an  hour  from  Fort  Collins  or  Boulder,  and  hour  and  a 
half  from  Denver  International  Airport,  is  one  of  the  most  beautiful  resort  towns  in  the  US,  with 
numerous  galleries,  shops,  and  restaurants,  along  with  majestic  nature  and  abundant  wildlife. 

The  shuttle  transportation  from  Denver  International  Airport  (DIA)  to  Estes  Park  (Stanley  Hotel) 
and  back  to  DIA  (the  shuttle  takes  about  lh  40  mins)  will  be  provided  by  Estes  Park  Shuttle. 
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Some  Details  about  Social  Program 


Trip  to  Rocky  Mountain  National  Park 

Wednesday  June  6,  2012,  8:00-14:00;  Picnic  lunch  included 

Rocky  Mountain  National  Park,  http://www.nps.gov/romo. 
established  as  NP  in  1915,  is  one  of  the  most  visited  national  parks  in 
the  US.  The  bus  tour  of  the  park  will  feature  a  guided  and  narrated 
drive  along  Trail  Ridge  Road  (the  highest  continuous  motorway  in  the 
US),  visits  to  Continental  Divide  and  the  start  of  the  mighty  Colorado 
River  (which  makes  the  Grand  Canyon  on  its  way  west),  to  the  “Land 
Above  the  Trees”  [reaching  12,183  feet  (3,713  m)  above  sea],  with  a 
lot  of  snow  (and  sun)  in  June  (see  http://www.rmnp.com),  a  visit  and 
longer  stay  at  the  Alluvial  Fan  http://www.rmnp.com/RMNP-Areas- 
HorseshoePark-AlluvialFan.html,  where  we  will  have  our  picnic 
lunch,  and  a  drive  (and  brief  stays)  through  lowland  meadows  (with 
herds  of  elk)  and  aspen  groves. 


Stanley  Hotel  Ghost  &  History  Tour 


Monday  17:30-19:00  tour  or  Tuesday  17:30-19:00  tour 

The  extremely  popular  spooky  sojourn  through  time  at  The  Stanley  Hotel, 
http ://stanlevho tel.com/tours/.  led  by  the  professional  and  experienced  tour  guides,  includes: 
1900’s  history  tour  featuring  the  beginnings  of  the  Hotel,  F.O.  &  Flora  Stanley  and  the  Stanley 
Steamer,  Stephen  King's  The  Shining  connection  to  the  Hotel  -  his  room  #217,  where  the 
creation  of  The  Shining  began,  ghost  stories  and  sightings,  the  Stanley  most  haunted  rooms  and 
places,  tours  through  the  underground  tunnel.  The  Ghost  &  History  Tour  lasts  approximately  90 
minutes. 


FEM2012  Closing  Reception 

Wednesday,  June  6,  2012,  18:00-21:00 

The  Elkhom  Lodge  &  Ranch,  Estes  Park  -  featuring  Chuckwagon  Dinner  and  Cowboy  Show 

Chuckwagon  Dinner  is  a  unique  Western  experience,  excellent  food  in  the  true  Old  West  setting, 
including  metallic  plates  etc.,  plus  some  cowboy  singers  and  entertainers.  The  Elkhom  Lodge  is 
the  oldest  continually  occupied  structure  in  Colorado. 
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About  Our  Financial  Sponsors 


ANSYS,  Inc. 


Today's  business  environment  is  rife  with  competitive  challenges,  customer  requirements  and 
financial  pressures.  These  factors  result  in  the  need  to  find  new  methods  for  engineering  more 
innovative  products  and  manufacturing  processes  —  while  minimizing  costs  and  time  to  market. 
Virtually  every  industry  now  recognizes  that  a  key  strategy  for  success  is  to  incorporate 
computer-based  engineering  simulation  early  in  the  development  process,  allowing  engineers  to 
refine  and  validate  designs  at  a  stage  where  the  cost  of  making  changes  is  minimal. 

The  ANSYS  simulation  software  portfolio  is  ideal  for  the  design  of  high-performance  electronic 
equipment  and  electromechanical  devices.  Companies  throughout  the  world  rely  on  ANSYS 
electromagnetic,  computational  fluid  dynamics  and  mechanical  analysis  solutions  when 
designing  mobile  communication,  internet  access,  broadband  networking  components  and 
systems,  integrated  circuits  (ICs)  and  printed  circuit  boards  (PCBs)  as  well  as  electromechanical 
systems  such  as  automotive  components  and  power  electronics  systems. 

At  ANSYS,  we  bring  clarity  and  insight  to  customers'  most  complex  design  challenges  through 
fast,  accurate  and  reliable  simulation.  Our  technology  enables  organizations  to  predict  with 
confidence  that  their  products  will  thrive  in  the  marketplace.  They  trust  our  software  to  help 
ensure  product  integrity  and  drive  business  success  through  innovation.  Founded  in  1970, 
ANSYS  employs  more  than  2,000  professionals,  and  many  of  them  are  expert  in  fields  such  as 
finite  element  analysis,  electronics  and  electromagnetics,  computational  fluid  dynamics  and 
design  optimization. 

Every  product  is  a  promise  to  live  up  to  and  surpass  expectations.  By  simulating  early  and  often 
with  ANSYS  software,  our  customers  become  faster,  more  cost-effective  and  more  innovative, 
realizing  their  own  product  promises. 


ANSYS 
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with 


For  more  information,  visit  www.ansys.com/HFSS 


Or  call  us  at  866  267  9724 


ANSYS  HFSS,  the  industry  standard  for  3D  full  wave  electromagnetic  field 
simulation,  is  now  accessible  through  an  easy  to  use  ECAD  based  user  in¬ 
terface.  This  new  and  fully  integrated  capability  will  allow  users  to  select, 
critical  structures,  nets,  transitions  in  the  layout  tool  and  easily  create 
models  to  be  solved  by  HFSS.  The  solution  setup  is  highly  automated,  easy 
to  use,  and  provides  the  high  fidelity  results  that  are  expected  with  HFSS. 
You  can  now  extract  highly  accurate  and  robust  S-parameters  model  with¬ 
out  ever  leaving  your  layout  environment.  This  same  capability  is  also 
available  from  within  the  Cadence  design  tools.  Ask  about  how  HFSS  for 
ECAD  can  enhance  your  SI  design  flow. 

U*0.'  r.  2  iil 

Realize  Your  Product  Promise™ 


About  Our  Financial  Sponsors 


CST 


CST  develops  and  markets  high  performance  software  for  the  simulation  of  electromagnetic 
fields  in  all  frequency  bands.  Its  products  allow  you  to  characterize,  design  and  optimize 
electromagnetic  devices,  all  before  going  into  the  lab  or  measurement  chamber.  This  can  help 
save  substantial  costs,  especially  for  new  or  cutting  edge  products,  reduce  design  risk,  and 
improve  overall  performance  and  profitability. 

Its  success  is  based  on  the  implementation  of  leading  edge  technology  in  a  user-friendly 
interface.  The  extensive  range  of  tools  integrated  in  CST  STUDIO  SUITE  enables  a  wide  range 
of  applications  to  be  analyzed  without  leaving  the  familiar  CST  design  environment.  This 
complete  technology  approach  includes  CST  MICROWAVE  STUDIO,  which  offers,  in  addition 
to  its  market  and  technology  leading  time  domain  solver,  FEM,  MoM,  and  TLM  solvers, 
enabling  unprecedented  simulation  reliability  through  cross  verification. 

CST  promotes  the  best-in-class  approach.  By  specializing  in  developing  3D  EM  software  and 
providing  straight-forward  links  with  other  best-in-class  vendors,  all  available  expertise  is  easily 
connected.  A  wide  range  of  import/export  filters  enable  the  easy  exchange  of  geometrical  data 
with  CAD  tools.  Furthermore,  imported  structures  can  be  modified  and  parameterized,  and  used 
for  optimization  and  design  studies.  Moreover  the  powerful  VBA  based  and  OLE  compatible 
macro  language  allows  direct  communication  with  programs  such  as  MATLAB. 

CST’s  customers  operate  in  industries  as  diverse  as  Telecommunications,  Defense,  Automotive, 
Electronics,  and  Medical  Equipment,  and  include  market  leaders  such  as  IBM,  Intel,  Mitsubishi, 
Samsung,  and  Siemens.  CST  markets  its  products  worldwide  through  a  network  of  distribution 
and  support  centers  which  also  provide  comprehensive  customer  support  and  training.  More 
information  is  available  at  http://www.cst.com. 
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Making  Sense  of  the  Real  World 
-  System  Level  EM  Simulation 


Simulation  of  cancer  treatment 
by  RF  thermoablation 


■  Components  don't  exist  in  electro¬ 
magnetic  isolation.  They  influence  their 
neighbors'  performance. They  are  affect¬ 
ed  by  the  enclosure  or  structure  around 
them.  They  are  susceptible  to  outside 
influences.  With  System  Assembly  and 
Modeling,  CST  STUDIO  SUITE  2012 
helps  optimize  component  as  well  as 
system  performance. 


If  you're  more  interested  in  filters,  cou¬ 
plers,  planar  and  multilayer  structures, 
we've  a  wide  variety  of  worked  applica¬ 
tion  examples  live  on  our  website  at 

www.cst.com/apps. 

Get  the  big  picture  of  what's  really  going 
on.  Ensure  your  product  and  components 
perform  in  the  toughest  of  environments. 


Involved  in  biomedical  applications?  You 
can  read  about  how  CST  technology  was 
used  to  simulate  biomedical  devices  at 

www.cst.com/biomed. 


Choose  CST  STUDIO  SUITE  2012  - 
complete  technology  for  3D  EM. 


CST 


CHANGING  THE  STANDARDS 


CST  of  America®,  Inc.  |  To  request  literature  (508)  665  4400  |  www.cst.com 


About  Our  Financial  Sponsors 


FEKO 

COMPREHENSIVE  ELECTROMAGNETIC  SOLUTIONS 

FEKO  (www.feko.info)  is  a  comprehensive  computational  electromagnetics  (CEM)  code  based 
on  accurate  “full  wave”  and  “asymptotic”  solutions.  Implementation  of  MLFMM  solver  along 
with  hybridization  of  MoM  with  asymptotic  high  frequency  techniques  (PO  and  UTD)  allows 
simulation  of  electrically  large  problems,  e.g.  antenna  placement.  Several  techniques  including 
Finite  Element  Method  (FEM)  are  available  for  dielectric  regions.  Examples  of  analysis  of 
electromagnetic  problems  include  EMC,  shielding,  coupling,  antenna  design,  antenna  placement 
analysis,  microstrip  antennas,  and  circuits,  striplines,  dielectric  media,  scattering  analysis  etc. 
The  implementation  is  comprehensive  and  has  been  extended  for  the  analysis  of  thin  dielectric 
sheets,  multiple  homogeneous  dielectric  bodies  and  planar  stratified  media.  Modules  are 
available  for  parametric  optimization,  time  domain  analysis  and  the  analysis  of  complex  cable 
harnesses  in  their  3D  environment.  FEKO  is  a  product  of  EM  Software  &  Systems  (SA)  Pty  Ltd. 

Contact: 

EM  Software  &  Systems  (USA)  Inc. 

144  Research  Drive 
Hampton  VA  23666 
Phone:  757-224  2035 
Fax:  757-282  5897 
Email:  feko@emssusa.com 
www.emssusa.com 
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Electrical  Size 


Complexity  of  Materials 


One  Product.  Multiple  Solvers. 

FEKO  includes  several  computational  methods,  each  optimised  for  different  problem  types.  Due  to  a  long 
history  of  hybridising  different  techniques,  FEKO  has  been  at  the  forefront  of  the  efficient  analysis  of 
complex,  low  and  high  frequency  problems.  The  Method  of  Moments,  Finite  Element  Method,  Multilevel 
Fast  Multipole  Method,  Uniform  Theory  of  Diffraction,  Physical  Optics  and  Ray-Launching  Geometrical 
Optics  are  all  available  in  the  standard  package. 

Additional  Applications:  Antenna  Design,  Antenna  Placement,  Waveguide,  RF  Components,  Microstrip 
Circuits,  EMC,  Cable  Coupling,  Radomes,  RCS,  Bio-EM. 


www.feko.info 


Global  sales  and  technical  support  network: 

Local  distributors  in  Europe,  North  America,  South  America,  Japan,  China,  South  Korea,  Singapore,  India,  Israel,  Taiwan,  South  Africa 


Detailed  Technical  and  Social  Program 


Sunday,  June  3,  2012 

19:00-21:00 

Welcome  reception  (Music  Room) 

Monday,  June  4,  2012 

7:00 

Breakfast  (MacGregor  Room) 

Billiard  Room 

Music  Room 

8:00 

Session  1:  Modeling  and  Design  of  Antennas  and  Arrays 

Session  2:  Discontinuous  Galerkin  Methods  (Q.  H.  Liu, 

Using  FEM  (R.  Kindt,  D.  Filipovic,  E.  Topsakal) 

W.  Cai) 

1.1.  Highly  Coupled  Bow-Tie  Array  Design  with  Passive 
Elements 

2.1.  On  the  Physics,  Accuracy,  and  Efficiency  of  the  Dual- 
Field  Domain-Decomposition  Method  for  Time-Domain 

Safa  Salman*,  Dimitris  Psychoudakis,  and  John  L.  Volakis 

Electromagnetic  Simulation 

Xiaolei  Li  and  Jian-Ming  Jin* 

8:20 

1.2.  Fully  Overlapping  Domain  Decomposition  in  Simulating 
Irregularly  Arranged  Arrays  in  Large  Media 

Tao  Peng*,  Kubilay  Sertel,  and  John  L.  Volakis 

2.2.  A  Full  Vectorial  Generalized  Discontinuous  Galerkin 
Beam  Propagation  Method  for  Optical  Waveguides 

Wei  Cai* 

8:40 

1.3.  Modeling  of  Antennas  and  Arrays  Using  Domain 
Decomposition  Method 

2.3.  Non-Conformal  and  Parallel  and  Interior  Penalty 
Discontinuous  Galerkin  Method  for  the  Time-Domain 

Kezhong  Zhao*  and  Nancy  Lambert 

Maxwell's  Equations 

Stylianos  Dosopoulos  and  Jin-Fa  Lee 

9:00 

1.4.  Theoretical  Analysis  of  Varactor-loaded  Half-width  Leaky- 
wave  Antenna  Arrays 

2.4.  Discontinuous  Galerkin  Methods  for  Ideal  MHD 

Equations 

Leo  C.  Kempel*,  B.  Shanker,  Edward  J.  Rothwell,  and  Prem 
Chahal 

Fengyan  Li* 

9:20 

1.5.  Non-conformal  Domain  Decomposition  Methods  for 
Modeling  Large  Finite  Antenna  Arrays 

Zhen  Peng,  Kheng-Hwee  Lim,  and  Jin-Fa  Lee 

2.5.  New  Spectral-Prism-Element  for  Layered  Structures 
implemented  in  an  efficient  implicit  DG-FETD  method 

Luis  E.  Tobon  and  Qing  H.  Liu* 

9:40 

Coffee  break  (Pinion  room) 

10:10 

1.6.  Investigation  of  the  Scattering  Behavior  of  Transponder 
Antennas  in  Case  of  Nonlinear  Termination 

2.6.  An  Efficient  Discontinuous  Galerkin  Finite-Element 
Time-Domain  Method  and  its  Application  to  Micro-Resistivity 

Thomas  Bauernfeind*,  Kurt  Preis,  Gergely  Koczka,  and  Oszkar 
Biro 

Imager 

Jiefu  Chen,  Carlos  Haramboure,  Lance  Pate,  Shawn 

Wallace,  Jay  Martin,  Medhat  Mickael,  and  Mac  Wisler 

10:30 

1.7.  On  the  Design  of  Finite  Planar  Ultrawideband  Modular 
Antenna  (PUMA)  Arrays  using  Domain  Decomposition  Methods 

Session  4:  FEM  and  CEM  Applications  in  Optics  and 
Nanophotonics  (E.  Simsek) 

Steve  S.  Holland,  Jack  T.  Logan,  Rick  W.  Kindt,  and  Marinos  N. 
Vouvakis* 

4.1.  Electromagnetic  Scattering  of  Finite  Size  Periodic  Large 
Scale  System:  Simulation  by  Domain  Decomposition 

Spectral  Element  Method 

Ma  Luo  and  Qing  Huo  Liu* 

10:50 

Session  3:  Time-Domain  FEM  and  Applications  (D.  White,  O. 
Biro,  T.  Rylander) 

4.2.  Efficient  alternatives  to  edge  vector  finite  elements  for 

EM  field  calculations  using  Hermite  interpolation  polynomial 

3.1.  Numerical  Simulation  of  HPM  Devices  with  ICEPIC 

basis  functions  on  a  nodal  grid  in  2D 

Andrew  D.  Greenwood* 

J.  D.  Albrecht,  C.  R.  Boucher,  C.  1.  Ahheng,  and  L.  R.  Ram- 
Mohan 

16 


Monday,  June  4,  2012 

Billiard  Room  Music  Room 


11:10  3.2.  Time  Domain  Method  for  Periodic  Nonlinear  Eddy  Current 
Problems 

Oszkar  Biro,  Gergely  Koczka,  and  Kurt  Preis 

4.3.  Parallel  Finite  Element  Electromagnetic  Code  Suite  for 
Accelerator  Modeling  and  Simulation 

Cho  Ng,  Lixin  Ge,  Kwok  Ko,  Kihwan  Lee,  Zenghai  Li,  Greg 
Schussman,  and  Liling  Xiao 

1 1 :30  3.3.  Hybrid  FETD/FDTD  Techniques  Based  on  the 
Discontinuous  Galerkin  Method 

Bao  Zhu,  Jiefu  Chen,  Wanxie  Zhong,  and  Qing  Huo  Liu 

4.4.  Optical,  Electromagnetics,  and  Thermal  Modeling  of 
Interaction  of  a  Focused  Beam  of  Light  with  Plasmonic 
Nanoparticles 

Eren  S.  Unlu  and  Kursat  Sendur* 

1 1 :50  Lunch  break  (MacGregor  Room) 


13:10  3.4.  Transient  Thermal  Analysis  using  a  Non-conformal  Domain  4.5.  Finite  Element  Analysis  of  Wideband  Nanostructures  for 
Decomposition  Approach  Photovoltaic  Applications 


Yang  Shao* ,  Zhen  Peng,  and  Jin-Fa  Lee 

Niccolo  Breda,  Giuseppe  Pelosi,  Stefano  Selleri,  and 

Ruggero  Taddei 

13:30  3.5.  An  Explicit  and  Unconditionally  Stable  Time-Domain  Finite- 
Element  Method  of  Linear  Complexity 

Qing  He,  Duo  Chen,  and  Dan  Jiao 

4.6.  Hybrid  Optical  Waveguides 

Ergun  Simsek*  and  Qing  H.  Liu 

13:50  3.6.  Conformal  PML  Modeling  in  DGTD  using  Continuous 
Material  Properties 

Jue  Wang*,  Zhen  Peng,  and  Jin-Fa  Lee 

Session  6:  Advances  in  Vector  Bases  for  CEM  (R. 

Graglia,  A.  Peterson,  D.  Wilton) 

6.1.  On  the  Construction  of  Well-conditioned  Hierarchical 
Bases  for  Tetrahedral  H(curl)-Conforming  Nedelec  Element 
Wei  Cai*  and  Jianguo  Xin 

14:10  Session  5:  Model  Order  Reduction  Techniques  (R.  Dyczij- 
Edlinger,  V.  de  la  Rubia) 

5.1.  Reduced-order  Electromagnetic  Modeling  in  the  Presence 
of  Uncertainty 

Juan  S.  Ochoa  and  Andreas  C.  Cangellaris* 

6.2.  High  Order  Finite  Elements  for  Computational  Physics: 

An  LLNL  Perspective 

Robert  N.  Rieben* 

14:30  5.2.  Efficient  Model  Order  Reduction  for  FEM  Analysis  of 
Waveguide  Structures  and  Resonators 

Fotyga  G.*,  Kulas  L.,  Nyka  K.,  and  Mrozowski  M. 

6.3.  FEM-MoM-Diakoptic  Analysis  of  Scatterers  with 
Anisotropic  Inhomogeneities  Using  Hierarchical  Vector 

Bases  on  Large  Curved  Elements 

Ana  B.  Manic*,  Dragan  1.  Olcan,  Milan  M.  Hie,  and  Branislav 
M.  Notaros 

14:50  5.3.  Model  Order  Reduction  for  Systems  with  Geometrical 
Parameters 

Kynthia  K.  Stavrakakis*,  Til man n  Wittig,  Wolfgang  Ackerman n, 
and  Thomas  Weiland 

6.4.  Rapid,  High-Order  Finite  Element  Modelling  with 

FEniCS  and  SUCEM:FEM 

A.J.  Otto*,  E.  Lezar,  N.  Marais,  R.G.  March  and,  and  D.B. 
Davidson 

15:10  Coffee  break  (Pinion  room) 


15:40  5.4.  Reduced-Order  Models  for  Geometric  Parameters 

Stefan  Burgard*,  Ortwin  Fade,  and  Romanus  Dyczij-Edlinger 

6.5.  Issues  Arising  in  Connection  with  Singular  Vector  Bases 
Roberto  D.  Graglia  and  Andrew  F.  Peterson* 

16:00  5.5.  Seamless  Reduced  Basis  Element  Method  for  Waveguide 
Simulation 

Yanlai  Chen*,  Jan  S.  Hesthaven,  and  Yvon  Maday 

6.6.  An  Orthogonal  Prism  Vector  Bases  Based  Quadratic 
Eigenvalue  Solver  of  Linear  Complexity  for  3-D 
Electromagnetic  Analysis 

Jongwon  Lee,  Duo  Chen,  and  Dan  Jiao 

16:20  5.6.  The  Reduced  Basis  Method  for  Saddle  Point  Problems 
Karen  Veroy*  and  Anna-Lena  Gerner 

6.7.  Preliminary  Research  on  the  Discontinuous  Enrichment 
Method  Based  Domain  Decomposition  Scheme  for  Solving 
the  Three-Dimensional  Vector  Curl-Curl  Equation 

Ming-Feng  Xue*  and  Jian-Ming  Jin 

16:40  5.7.  Reliable  Fast  Frequency  Sweep  of  Microwave  Circuits  via 
the  Reduced-Basis  Method 

Valentin  de  la  Rubia 

6.8.  An  Iterative  Domain  Decomposition  Method  for  Solving 
Wave  Propagation  Problems 

Gergely  Koczka,  Kurt  Preis*,  Thomas  Bauernfeind,  and 
Oszkar  Biro 

17:00  5.8.  Reduced  Order  Modeling  of  High  Magnetic  Field  Magnets 

C.  Daversin,  C.  Prudhomme,  C.  Trophime*,  and  S.  Veys 

6.9.  Evaluating  Singular  and  Near-Singular  Integrals  on 
Curvilinear  Elements 

Donald  R.  Wilton* 

17:30-  19:00  Stanley  Hotel  Ghost  &  History  Tour 

19:30-  22:30  Conference  Banquet  (MacGregor  Room)  (the  room  opens  at  19:00) 
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7:00  Breakfast  (MacGregor  Room) 

Billiard  Room  Music  Room 


Session  7:  FETD  Modeling  of  Complex  Media  and 

Session  8  Adaptive  FEM,  Higher  Order  Bases,  and 

8:00  Structures  (R.  Lee,  F.  Teixeira) 

7.1.  FETD/FDTD  with  Automatic  Hybrid  Mesh  Generation 

Shu  min  Wang* 

Advanced  FEM  Formulations  (J.  Webb,  M.  Gavrilovic,  M. 
Salazar  Palma,  L.  G.  Castillo) 

8.1.  Hp-adaption  for  Computing  Scattering  Parameters  on 
Tetrahedral  Meshes 

J.  P.  Webb 

8:20  7.2.  Study  of  Numerical  Dispersion  for  Quadrilateral  and 

Triangular  Elements  Used  in  the  Mixed  Finite  Element  Method 
Luis  E.  Tobon  and  Qing  H.  Liu* 

8.2.  Efficient  RCS  Prediction  of  Jet  Engine  Air  Intakes  Using 
hp-FEM  and  Reduced  Order  Modeling 

Adam  Zdunek*  and  Waldemar  Rachowicz 

8:40  7.3.  Probabilistic  Finite-Difference  Time-Domain  Method 
Through  Stochastic  Electromagnetic  Macro-Models 

Ata  Zadehgol*,  Andreas  C.  Cangellaris 

8.3.  Rules  for  Adoption  of  Expansion  and  Integration  Orders 
in  FEM  Analysis  Using  Higher  Order  Hierarchical  Bases  on 
Generalized  Hexahedral  Elements 

Nada  J.  Sekeljic*,  Slobodan  V.  Savic,  Milan  M.  Ilic.  and 

B  rani  Slav  M.  Notaros 

9:00  7.4.  Advancing  Ultrahigh  Field  Human  MRI  using  FDTD  Method  8.4.  Auxiliary  Space  Preconditioning  for  Hierarchical 
TamerS.  Ibrahim*  Elements 


A.  Aghabarati,  J.  P.  Webb 

9:20  7.5.  The  Use  of  Multiple  Time  Steps  in  the  FDTD  Method 

Ruben  Ortega  and  Robert  Lee* 

8.5.  Chain  Homotopy  Operators  and  Whitney  form-based 
Finite  Element  Analysis 

P.  Robert  Kotiuga* 

9:40  Coffee  break  (Pinion  room) 


10:10  Session  9:  Advances  in  Hybrid  Methods  and  Multiphysics 
Problems  (B.  Shanker,  L.  Kempel) 

9.2.  Thermal-Aware  DC  IR-Drop  Co-Analysis  using  a  Non- 
conformal  Domain  Decomposition  Approach 

Yang  Shao*,  Zhen  Peng,  Jin-Fa  Lee 

8.6.  Comparative  Performance  of  Nodal-Based  Versus 
Edge-Based  Finite  Element  Formulations 

Ruben  Otin*,  Santiago  Badia,  and  Luis  E.  Garcia-Castillo 

10:30  9.3.  Development  of  Time  Domain  Discontinuous  Galerkin- 
Vector  Generalized  Finite  Element  Method 

O.  Tuncer,  B.  Shanker*,  L.  C.  Kempel 

8.8.  A  Three-Dimensional  Self-Adaptive  hp  Finite  Element 
Method  for  the  Characterization  of  Waveguide 

Discontinuities 

Ignacio  Gomez-Revuelto,  Luis  E.  Garcia-Castillo,  Sergio 
Llorente-Romano,  and  David  Pardo 

10:50  9.4.  Multi-solver-based  Generalized  Impedance  Boundary 
Condition  for  Complicated  EM  Simulation 

Shiquan  He,  Jun  Hu*,  Zaiping  Nie,  Lijun  Jiang,  W.  C.  Chew 

Session  10:  FEM  Modeling  and  Applications  of 
Metamaterials  and  Periodic  Media  (K.  Sertel,  J.  Volakis) 

10.1.  Finite  Element  Modeling  of  Magnetic  Metamaterials 
and  Associated  Challenges  for  Their  Experimental 
Characterization 

Nil  Apaydin* ,  Kubilay  Sertel,  and  John  L.  Volakis 

11:10  9.5.  Analysis  of  Electromechanical  Devices  Using  the  Dual- 
Primal  Finite  Element  Tearing  and  Interconnecting  Method 
Incorporated  with  the  LU  Recombination  Method 

Wang  Yao*,  Jian-Ming  Jin 

10.2.  A  Method  for  Rapid  Evaluation  of  Space-time 
Convolutions  with  the  Time  Domain  Periodic  Green's 

Function 

D.  Dault  and  B.  Shanker* 

1 1 :30  9.6.  Alternative  Formulations  for  Hybrid  Electromagnetic-Circuit  10.3.  Finite  Element  Analysis  of  Photon  Density  of  States  for 
Simulators  Two-Dimensional  Photonic  Crystals  with  Omnidirectional 


Vivek  Subramanian,  AH  E.  Yilmaz* 

Light  Propagation 

R.  F.  Jao  and  M.  C.  Lin 

1 1 :50  Lunch  break  (MacGregor  Room) 


13:10  9.7.  Symmetric  and  Nonsymmetric  FEM-MoM  Techniques  Using  10.4.  A  Hybrid  Approach  for  Characterizing  Large 

Higher  Order  Hierarchical  Vector  Basis  Functions  and  Curved  Anisotropic  and  Inhomogeneous  Metamaterial  Structures 


Parametric  Elements 

Ana  Manic*,  Milan  Ilic,  B  rani  Slav  Notaros 

Justin  G.  Pollock  and  Ashwin  K.  Iyer* 
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Music  Room 


13:30  Session  11:  Domain-Decomposition  Methods  (J.-F.  Lee,  Z.  Session  12:  Advanced  FEM/MoM  Modeling,  Design,  and 
Peng)  Optimization  (C.J.  Reddy,  U.  Jakobus,  J.  Zapata,  J.  Gil, 

11.1.  Analysis  of  Three-Dimensional  Array  Structures  Using  D.  Jiao,  A.  Cangellaris) 

Nonconformal  and  Cement  FETI-DP  Methods  12.1.  Overview  of  the  Hybrid  Finite  Element  Method/Method 

Ming-Feng  Xue*  and  Jian-Ming  Jin  of  Moments  Capability  in  FEKO 

_ Marianne  Bingle*,  Ulrich  Jakobus,  and  Johann  J.  van  Tonder 

13:50  11.2.  A  Domain-Decomposition-Based  Preconditioner  of  FE-BI-  12.2.  Application  of  Hybrid  FEM/MoM  Technique  for 
MLFMAfor  Computing  EM  Scattering  by  Large  Inhomogeneous  Handheld  Mobile  Devices 
3D  Targets  Rensheng  Sun  and  C.  J.  Reddy 

_ Ming-Lin  Yang,  Xin-Qing  Sheng,  and  Cheng-Dan  Huang _ 

14:10  11.3.  On  the  Computational  Complexity  of  Domain  Wise  12.3.  Application  of  Hybrid  FEM/MoM  Technique  for  Cavity 

Multilevel  Multifrontal  Method  Backed  Apertures 

_ Peng  Liu*  and  Chao-Fu  Wang _ Gopinath  Gampala,  Rensheng  Sun  and  C.  J.  Reddy _ 

14:30  11.4.  Harmonic  Balance  Domain  Decomposition  Finite  Element  12.4.  Eigenproblem  Approach  for  Analysis  and  Optimization 

for  Nonlinear  Passive  Microwave  Devices  Analysis  of  Microwave  Filters  Using  Finite  Element  Method 

_ Laurent  Ntibarikure,  Giuseppe  Pelosi,  and  Stefano  Selleri _ Adam  Lamecki,  Michal  Mrozowski _ 

14:50  11.5.  Application  of  Tree-Cotree  Splitting  to  the  Dual-Primal  12.5.  Does  Sensitivity  Analysis  for  Fast  Frequency  Sweeps 

Finite  Element  Tearing  and  Interconnecting  Method  for  Solving  Converge? 

Low-Frequency  Breakdown  Problems  Andreas  Kohler,  Ortwin  Fade,  Romanus  Dyczij-Edlinger* 

_ Wang  Yao*  and  Jian-Ming  Jin _ 

15:10  Coffee  break  (Pinion  room) 

15:40  11.6.  Edge  Elements  and  DPM  with  PBSV  to  Accelerate  Solving  12.6.  Benchmark  of  Conformal  Finite-Difference  Method 
Scattering  Problems  of  Electrically  Large  Objects  against  Finite  Element  Method 

Lianyou  Sun*,  Lishan  Xue,  and  Wei  Hong_ M.  C.  Lin 


16:00  11.7.  Enhancement  of  Second-harmonic  Generation  in  an  Air-  12.7.  Layout-Integrated  Electromagnetic  Interconnect 
bridge  Photonic  Crystal  Slab:  Simulation  by  Spectral  Element  Characterization  and  Simulation 
Method  Dennis  Nagle,  Jilin  Tan,  Jian-Ming  Jin 


_ Ma  Luo  and  Qing  Huo  Liu* _ 

16:20  11.8.  A  Non-overlapping  and  Non-conformal  Domain  12.8.  An  Extraction-Free  Circuit  Simulator  of  Linear 

Decomposition  Method  with  Optimized  Second  Order  Complexity  Guided  by  Electromagnetics-Based  First 

Transmission  Conditions  for  Time-Harmonic  Maxwell  Equations  Principles 

in  R  Qing  He,  Duo  Chen,  and  Dan  Jiao 


Zhen  Peng  and  Jin-Fa  Lee 


16:40  11.9.  The  Quest  for  Robust  Domain  Decomposition  Methods  in  12.9.  A  model  reduction  method  for  multiple  particle 
CEM:  Are  We  There  Yet?  electromagnetic  scattering  in  three  dimensions 

_ Marinos  N.  Vouvakis*  and  Georgios  N.  Paraschos _ M.  Ganesh*,  J.  Hesthaven,  and  B.  Stamm _ 

17:00  Technical  program  ends  for  the  day 

17:00-  18:30  FEM2012  Scientific  Committee  and  Organizing  Committee  Meeting  (Music  Room) 

17:30-  19:00  Stanley  Hotel  Ghost  &  History  Tour 


Wednesday,  June  6,  2012 

7:00-  8:00 

Breakfast  (MacGregor  Room) 

8:00-  14:00 

Trip  to  Rocky  Mountain  National  Park  (picnic  lunch) 

Open  Forum  Discussions  (Music  Room) 

15:00  Open  forum  discussion  1:  Open  Forum  Discussion  on  FEM  and  CEM  Bases,  Elements,  and  Formulations 


_ (A.  Peterson,  D.  Wilton,  J.  Webb,  K.  Sertel) _ 

15:45  Open  forum  discussion  2:  Open  Forum  Discussion  on  DDM,  MOR,  and  Efficient  FEM  and  Hybrid  Solutions 
_ (J.-F.  Lee,  Zhen  Peng,  R.  Dyczij-Edlinger,  V.  de  la  Rubia) _ 

16:30 _ Coffee  break  (Pinion  Room) _ 

16:50  Open  forum  discussion  3:  Open  Forum  Discussion  on  FEM-Based  Modeling,  Design,  and  Applications 
_ (L.  Kempel,  B.  Shanker,  R.  Kindt,  C.J.  Reddy,  K.  Zhao,  T.  Euler) _ 

18:00-  21:00_ Closing  Reception  (The  Elkhorn  Lodge  &  Ranch  -  Chuckwagon  Dinner  and  Cowboy  Show) 
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Billiard  Room 

Session  1:  Modeling  and  Design  of  Antennas  and  Arrays  Using  FEM 
(R.  Kindt,  D.  Filipovic,  E.  Topsakal) 

1.1.  Highly  Coupled  Bow-Tie  Array  Design  with  Passive  Elements 

Safa  Salman*,  Dimitris  Psychoudakis,  and  John  L.  Volakis 

8:20 

1.2.  Fully  Overlapping  Domain  Decomposition  in  Simulating  Irregularly 

Arranged  Arrays  in  Large  Media 

Tao  Peng*,  Kubilay  Sertel,  and  John  L.  Volakis 

8:40 

1.3.  Modeling  of  Antennas  and  Arrays  Using  Domain  Decomposition  Method 
Kezhong  Zhao*  and  Nancy  Lambert 

9:00 

1.4.  Theoretical  Analysis  of  Varactor-loaded  Half-width  Leaky-wave  Antenna 

Arrays 

Leo  C.  Kempel*,  B.  Shanker,  Edward  J.  Rothwell,  and  Prem  Chahal 

9:20 

1.5.  Non-conformal  Domain  Decomposition  Methods  for  Modeling  Large  Finite 
Antenna  Arrays 

Zhen  Peng,  Kheng-Hwee  Lim,  and  Jin-Fa  Lee 

9:40 

Coffee  break  (Pinion  room) 

10:10 

1.6.  Investigation  of  the  Scattering  Behavior  of  Transponder  Antennas  in  Case 
of  Nonlinear  Termination 

Thomas  Bauernfeind*,  Kurt  Preis,  Gergely  Koczka,  and  Oszkar  Biro 

10:30 

1.7.  On  the  Design  of  Finite  Planar  Ultrawideband  Modular  Antenna  (PUMA) 

Arrays  using  Domain  Decomposition  Methods 

Steve  S.  Holland,  Jack  T.  Logan,  Rick  W.  Kindt,  and  Marinos  N.  Vouvakis* 
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Highly  Coupled  Bow-Tie  Array  Design  with  Passive  Elements 

Safa  Salman*,  Dimitris  Psychoudakis,  and  John  L.  Volakis 
ElectroScience  Laboratory,  Dept,  of  Electrical  and  Computer  Engineering,  The  Ohio  State 
University,  Columbus  OH  43212,  USA  salman.9@osu.edu,  dpsycho@ece.osu.edu, 

volakis@ece.osu.edu 


Increasing  demand  for  high  speed  wireless  communication  systems  implies  a  need  for  broadband  planar 
antenna  arrays.  Broadband  conformal  antennas  become  preferable  for  both  transportation  and  installation 
on  platforms,  due  to  several  advantages  including  their  light  weight  and  low-cost  production. 

In  this  work,  the  concept  of  a  current  sheet  array  (CSA)  is  expanded  to  design  a  broadband  bowtie¬ 
shaped  antenna  array  using  HFSS.  A  bowtie-shaped  antenna  element  was  chosen  as  the  element  itself  has 
good  bandwidth  and  radiation  characteristics.  The  antenna  design  parameters  were  chosen  for  operating 
over  the  VHF  and  UHF  frequency  bands.  As  the  overlapping  array  elements  illustrated  in  the  figure 
suggest,  there  is  a  thick  substrate  layer  between  the  overlapping  elements  which  has  a  dielectric  constant 
of  sr=10.2.  A  key  concept  of  the  design  was  the  introduction  of  strong  mutual  coupling  between  the 
elements,  using  FEM  (via  HFSS)  to  optimize  the  spacing  between  the  elements.  The  present  approach  has 
the  potential  to  achieve  a  7:1  bandwidth  over  the  VHF  and  UHF  frequency  bands  with  wide-angle  beam. 

A  small  3x3  finite  array  was  simulated  using  HFSS  and  then  fabricated  and  experimentally  tested  to 
verify  that  the  broadband  characteristics  can  still  be  achieved  using  a  small  aperture  of  overlapped  bowtie 
dipoles.  The  physical  realization  of  the  proposed  antenna  array  is  also  shown  in  the  figure  below.  Both 
top  and  bottom  views  are  illustrated.  The  reflection  coefficient  of  the  finite  array  is  also  shown  in  the 
same  figure.  For  testing,  only  the  center  element  was  excited  while  the  others  were  terminated  with  75  Q. 
This  array  structure  was  backed  with  a  finite  ground  plane  of  the  same  dimensions  as  the  dielectric  slab 
and  the  array  both  for  FEM  solutions  and  measurements.  At  the  lowest  frequency  the  total  array  thickness 
was  X0/4.  As  the  comparison  between  simulated  and  measured  data  shows,  the  simulated  and  measured 
data  follow  each  other  well  however  more  fabrication  care  needed  to  insure  better  agreement. 

Future  work  may  involve  the  use  of  flat  antenna  array  realizations  for  more  than  7:1  operating  bandwidth 
design  target. 


Left:  Unit  cell  of  the  overlapping  bowtie  antenna  array  over  a  dielectric  slab  with  infinite  ground  plane.  Middle:  Center 
element  reflection  coefficient  with  all  other  elements  terminated  at  75  Q.  Right:  Physical  layout  of  the  printed  element  layer  of 
the  array  top  (left)  and  bottom  (right)  views. 
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A  fully  overlapping  domain  decomposition  (DDM)  algorithm  is  proposed  to  simulate  small  features 
within  large  domains.  A  unique  aspect  of  this  approach  is  that  it  does  not  require  a  common  interface 
between  separate  subdomains  to  glue  them  together.  Instead,  it  employs  a  coarse  background  domain  with 
several  independently  meshed  subdomains  to  model  embedded  details.  Much  like  other  domain 
decomposition  schemes,  this  approach  shares  similarities  to  Multigrid  methods,  particularly  with  adaptive 
refinement  Multigrid.  However,  most  Multigrid  methods  either  treat  the  coarse  domain  as  a 
preconditioning  tool  (V-cycle  Multigrid)  or  formulate  the  problem  as  a  local  h  refinement  scheme  (fast 
adaptive  composite  Multigrid).  The  coarse  background  meshes  in  these  setups  must  invariably  resolve  the 
small  features  to  obtain  correct  error  updates.  In  contrast,  our  approach  does  not  need  to  resolve  the  small 
features  within  the  coarse  background  domain.  Most  importantly,  the  introduced  multiple  detail  domains 
do  not  form  conforming  common  interfaces  as  in  traditional  non-overlapping  DDMs.  Instead, 
communication  between  the  subdomains  (modeling  the  details)  and  the  coarse  background  mesh  is  done 
by  treating  the  subdomain  fields  as  sources  to  the  background  domain.  These  flexibilities  enable  a 
modeling  ease  and  much  enhanced  convenience.  In  this  paper,  we  employ  the  proposed  non-overlapping 
DDM  to  model  arrays  that  may  include  non-periodic  elements. 

Modeling  highly  repetitive  structures,  as  is  the  case  with  arrays,  is  a  well  studied  process.  It  has  been 
considered  in  the  context  of  finite  element  (FEM)  and  integral  methods.  In  the  case  of  FEM -DDM  with 
sub-structuring,  multiple  subdomains  are  used  to  model  the  large  problem,  and  the  unknowns  are 
typically  those  at  the  interface  of  the  smaller  subdomains,  implying  significant  reduction  in  the  number  of 
unknowns.  To  glue  all  subdomains,  bi-directional  Robin  Transmission  Conditions  are  typically  used  on 
the  geometrically  conformal  domain  boundaries.  Typically,  each  subdomain  must  be  geometrically 
connected  to  all  its  neighbors.  Integral  methods,  that  may  also  incorporate  the  multilevel  FMM 
(MLFMM),  rely  on  the  same  divide-and-conquer  philosophy,  but  can  be  cumbersome  when 
inhomogeneous  media  are  involved.  In  contrast  to  the  above,  the  fully  overlapping  DDM  offers  a 
relatively  carefree  solution  due  to  its  freedom  and  flexibilities  in  modeling  antenna  small  features  and  in 
configuring  irregular  arrays.  Importantly,  trial  and  error  and  in  situ  optimizations  are  frequently 
indispensable  in  designing  antennas.  As  the  detail  subdomains  are  meshed  independently,  remeshing  of 
the  subdomains  during  optimizations  can  be  carried  out  with  ease.  Therefore,  the  proposed  non¬ 
overlapping  DDM  is  very  attractive  for  design  optimization.  In  this  paper,  the  fully  overlapping  DDM 
would  be  applied  to  multiple  irregularly  arranged  arrays  embedded  in  inhomogeneous  media.  This  is  done 
in  the  context  of  the  finite  element  method. 
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Computational  electromagnetics  plays  a  crucial  role  in  accurate  predictions  of  design  performance  prior 
to  construction.  Accurate  simulation  of  electrically  large  and  geometrically  complicated  electromagnetic 
(EM)  problems  is  of  vital  importance  in  many  areas  of  electrical  engineering,  but  also  is  a  very 
challenging  task.  The  scope  and  application  of  traditional  approaches  such  as  the  finite  element  method 
(FEM),  boundary  element  method  (BEM)  and  finite  difference  method  (FDM)  are  limited  to  moderate 
electrical  size  and  reduced  complexity.  These  limitations  stem  from  the  vast  computational  resources 
required  by  these  numerical  methods. 

The  domain  decomposition  method  (DDM)  has  emerged  as  a  powerful  and  attractive  technique  for  the 
analysis  of  large-scale  electromagnetic  problems  due  to  its  inherent  parallelism  and  its  beauty  as  an 
efficient  and  effective  preconditioner.  DDM  is  based  on  a  divide-and-conquer  philosophy.  Instead  of 
tackling  a  large  and  complex  problem  directly  as  a  whole,  the  original  problem  is  partitioned  into  smaller, 
possibly  repetitive,  and  easier  to  solve  sub-domains.  In  this  paper  DDM  is  used  as  an  effective  FEM 
preconditioner  where  a  higher  order  Robin’s  transmission  condition  (TC)  is  devised  to  enforce  the 
continuity  of  electromagnetic  fields  between  adjacent  sub-domains  and  at  the  same  time  accelerate  the 
convergence  of  the  iterative  process.  DDM  is  also  employed  to  formulate  a  novel  hybridization  of  FEM 
with  the  boundary  element  method  that  exact  treatment  of  the  radiation  condition  can  be  realized. 
Through  this  DDM-based  hybrid  finite  element-boundary  integral  (FE-BI)  method,  we  can  allow  FEM- 
domains  to  be  disconnected;  the  coupling  between  disjoint  domains  is  naturally  taken  into  account  via  the 
Green’s  function.  The  advantages  of  DDM-based  FE-BI  compared  to  traditional  FE-BI  include 
modularity  of  FEM  and  BI  domains  in  terms  of  mesh  and  selection  of  basis  functions.  This  “non- 
conformal”  ability  significantly  simplifies  the  integration  of  existing  state-of-art  FEM  and  BEM 
solvers.  The  continuity  enforcement  through  Robin’s  TC  naturally  renders  this  FE-BI  formulation  free  of 
internal  resonance  issues.  Since  domains  are  allowed  to  be  disjoint,  if  one  or  more  sub-domains  are  purely 
metallic  or  highly  conducting,  the  present  DDM  can  allow  integral  equation  methods  to  be  applied  to 
these  sub-domains  directly  so  as  to  reduce  memory  consumption.  Several  real-life  antennas  and  finite 
antenna  arrays  will  be  analyzed  to  demonstrate  the  power  of  this  method.  In  particular  we  will  show  that 
DDM  enables  scientists  and  engineers  to  solve  much  larger  problems  via  distributed  computing. 
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A  subject  of  investigation  over  the  past  few  years  involves  methods  to  control  the  radiation  properties  of 
half- width  leaky- wave  (HWLW)  antennas  using  reactive  tuning  elements.  HWLW  antennas  exhibit 
dispersion  in  both  the  voltage  standing- wave  ratio  (VSWR)  and  radiation  pattern.  The  impact  of  the 
former  can  be  mitigated  with  an  adaptive  matching  network  that  incorporates  suitable  reactive  elements. 
The  latter  can  also  be  mitigated,  over  a  more  limited  bandwidth,  through  the  use  of  reactive  loads  and  in 
doing  so,  the  requirements  for  a  VSWR  matching  network  are  reduced.  The  effect  of  pattern  dispersion  is 
that  the  main  beam  steers  from  near  broadside  to  near  end-fire  as  the  frequency  is  chirped.  Hence,  the 
antenna  “looks”  at  a  different  portion  of  the  scan  space  at  different  frequencies.  For  some  applications,  it 
is  desirable  that  the  antenna  has  a  stable  pointing 
direction  as  the  frequency  is  chirped. 


In  a  recent  presentation  (L.  Kempel,  E.  Rothwell,  B. 

Shanker,  and  P.  Chahal,  “Theoretical  analysis  of  a 
varactor-loaded  half-width  leaky-wave  antenna,”  XXX 
URSI  General  Assembly  and  Scientific  Symposium , 

Istanbul,  TURKEY,  Aug.  2011)  results  for  a  single 
element  are  presented  demonstrating  pattern  stability 
over  nearly  1  GHz  in  C-band.  Included  in  the  paper  is 
an  empirically  determined  second-order  dispersive 
compensation  load  equation: 

Cpf(  fam)  =  0.0032  -  0.0668  +  0.3228 

In  this  paper,  results  for  such  elements  placed  in  arrays 
will  be  presented.  Of  particular  interest  is  whether  the 
load  dispersion  requirements  change  due  to  mutual 
coupling.  The  rationale  for  this  line  of  inquiry  lies  in 
the  propagation  properties  of  the  leaky-wave  mode.  Since  it  has  both  a  longitudinal  and  in-plane 
orthogonal  components,  the  latter  may  provide  a  significant  coupling  mechanism  that  could  alter  the 
needed  reactive  compensation. 

The  analysis  will  be  performed  using  the  finite  element-boundary  integral  (FE-BI)  method.  It  is  well 
suited  for  this  investigation  since  it  permits  representation  of  dielectric  as  well  as  metallic  materials  in 
addition  to  reactive  lumped  (though  dispersive)  loads.  In  addition,  since  simulation  at  many  frequencies  is 
required,  the  memory  and  operation  count  efficiency  of  the  FE-BI  method  results  in  fairly  small 
calculation  times. 


Angle  [deg] 


Figure  1 .  Radiation  pattern  of  the 
varactor-loaded  HWLW  antenna  at  three 
frequencies. 
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Large  finite  antenna  arrays  are  commonly  used  in  wireless  communication  systems  and  radars  to 
transmit  and  receive  signals  through  spaces.  These  arrays  often  involve  complicated  and  multi¬ 
scale  geometrical  features  and  they  are  usually  electrically  large.  Thus,  full  wave  analysis  of  such 
problems  by  traditional  methods  demands  prohibitive  computer  resources. 

We  proposed  a  non-conformal  domain  decomposition  method  (NC-DDM)  for  solving  large  antenna 
arrays  mounted  on  platform,  see  Fig.l.  To  further  improve  the  convergence  in  the  DDM  iterations, 
an  optimal  2nd  order  transmission  condition  is  introduced  to  enforce  field  continuities  across  domain 
interfaces.  Moreover,  many  antenna  arrays  are  also  electrically  large.  Consequently,  the  use  of 
absorbing  boundary  condition  may  not  be  adequate  as  an  accurate  mesh  truncation  method.  Herein, 
we  combine  directly  the  finite  element  domain  decomposition  method  with  a  generalized  combined 
field  integral  equation  and  form  automatically  the  hybrid  finite  element  and  boundary  integral 
(FEBI)  method.  The  use  of  the  boundary  integral  method,  arguably,  offers  the  best  accuracy  for 
modeling  unbounded  electromagnetic  radiation  and  scattering  problems,  albeit  at  the  increases 
of  memory  and  CPU  times.  Furthermore,  the  finite  element  tearing  and  interconnecting  (FETI) 
method  is  employed  to  take  advantage  of  the  repetitions  to  drastically  reduce  the  computational 
resources. 


Figure  1:  An  interlaced  dual  periodic  array  mounted  on  an  aircraft 
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Introduction:  In  passive  backscattering  applications  like  UHF-RFID  (ultra  high  frequency-radio  frequency 
identification)  the  communication  between  the  interrogator  unit  (reader)  and  the  transponder  is 
established  by  means  of  modulating  the  radar  cross  section  of  the  transponder.  In  general,  for  UHF-RFID 
applications,  this  is  done  by  switching  the  analog  input  impedance  between  two  states  in  phase  with  the 
data  stream  to  be  transmitted.  Unfortunately  the  analog  input  impedance  is  not  constant,  indeed  it  has  a 
strong  nonlinear  behavior.  This  behavior  is  mainly  caused  by  a  voltage  limiter  at  the  transponder  IC’s 
frontend.  Hence,  the  transponder’s  input  impedance  is  a  function  of  the  induced  antenna  voltage.  In  the 
present  work  a  method  based  on  an  analytical  approach  (R.  C.  Hansen,  “Relationships  Between  Antennas 
as  Scatterers  and  as  Radiators,”  Proceedings  of  the  IEEE,  Vol.  77,  No.  5,  1989,  pp.  659-662)  is  used  to 
investigate  the  scattered  field  of  the  transponder.  The  analyses  are  carried  out  by  the  method  of  finite 
elements. 


Scattering  investigations:  In  general  the  scattered  field  is  investigated  by  means  of  a  full-wave  finite 
element  simulation  of  the  entire  communication  channel  including  the  antenna  of  the  reader,  the  air 
volume  between  the  reader  and  the  transponder  and  the  transponder  antenna.  Hence,  the  system  of 
equations  can  be  quite  large.  To  reduce  the  computational  effort  the  following  relationship  introduced  by 
Hansen  could  be  used: 


^scoA^l}  ^scaA^a}  .  F  Knt 

*  *  7 


(1) 


In  (1)  E scat(ZL)  is  the  scattered  field  of  the  transponder  if  the  transponder  antenna  is  connected  to  the  load 
impedance  ZL,  whereas  E scat(Za  )  is  the  scattered  field  if  the  transponder  IC  is  conjugate  complex  matched 
to  the  transponder  antenna’s  input  impedance  Za.  Eant  is  the  field  which  is  radiated  by  the  transponder 
antenna  if  the  antenna  is  driven  with  the  current  A,  T*  is  the  reflection  coefficient  and  Im  is  the  current 
which  flows  in  the  antenna  in  the  case  of  the  conjugate  complex  matching.  Hence  the  influence  of  the 
nonlinear  transponder  can  be  analyzed  by  simply  scaling  the  solution  of  the  driven  antenna  with  the 
reflection  coefficient  and  superimposing  the  solution  of  the  scattered  field  in  case  of  the  conjugate 
complex  matching.  Doing  so,  the  full-wave  finite  element  analyses  can  be  reduced  to  a  minimum  set.  If 
the  non-linear  voltage  dependent  behavior  of  the  transponder  IC  is  known,  e.g.  from  measurement,  one  is 
able  to  analyze  the  robustness  of  the  transponder  antenna  design  to  changes  in  the  IC  impedance. 
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The  planar  ultra- wideband  modular  array  (PUMA)  (S.  S.  Holland  et  ah,  IEEE 
Trans.  Antennas  Propag.,  Vol.  60,  p.  130,  2012)  is  a  new  class  of  wideband 
arrays  capable  of  5:1  bandwidth  from  a  totally-planar,  low-profile,  modular 
topology  using  a  feeding  technique  that  permits  direct  connection  to  standard 
unbalanced  50D  RF  interfaces.  Extensive  studies  on  infinite  array  unit  cells 
have  resulted  into  VSWR  levels  well  below  2,  but  very  work  little  has  been 
done  towards  understanding  the  impedance  and  far-field  radiation  properties 
under  realistic  finite  array  scenarios  (truncation  effects). 

This  work  will  investigate  the  radiation  properties  of  finite  PUMA  arrays  using 
a  domain  decomposition  finite  element  method  (DD-FEM)  (G.  N.  Paraschos  et 
al.  IEEE  Antennas  Prop.  Symp  2009,  DOI:  10. 1109/APS. 2009. 5171582).  The 
proposed  DD-FEM  is  based  on  the  finite  tearing  and  interconnecting  (FETI) 
method  with  two  Lagrange  multiplier  sets  per  interface,  effectively  enabling 
the  efficient  modeling  of  repetitive  structures  such  as  phased  arrays,  by  only 
storing  memory  for  one  array  element,  and  several  other  building  blocks  that 
represent  the  surrounding  array  region.  Even  with  this  DD-FEM  analysis,  a 
full  finite  array  simulation  is  computationally  demanding  and  not  appropriate 
for  design.  Alternatively,  a  faster  intermediate  modeling  step  is  introduced 
that  involves  solving  for  an  finite  x  infinite  array,  and  is  used  as  the  first-pass 
design.  This  is  achieved  by  combining  the  FETI-2A  with  periodic  boundary 
conditions  (PBCs)  (W.  Wang  et  al.  IEEE  Trans.  Antennas  Prop.  Vol.  59,  p. 
4142,  2011).  Once  dominant  trends  on  the  array’s  element  design  and  size  are 
established,  the  full  FETI-2A  DD-FEM  method  is  used  to  accurately  predict 
the  finite  PUMA  performance. 

This  talk  will  begin  by  introducing  the  basic  topology  of  the  PUMA  array 
highlighting  the  modeling  and  design  challenges  arising  from  the  multi-scale 
nature  of  the  geometry.  The  rational  and  the  basic  principles  behind  the  DD- 
FEM  will  be  presented.  A  baseline  test  case  of  a  3:1  bandwidth  16x16  PUMA 
will  be  used  to  present  results  for  the  embedded  element  impedance,  radiation 
pattern,  active  reflection  coefficient  and  patterns  under  various  exception  con¬ 
figurations.  The  results  will  be  compared  with  measurements,  and  trends  on 
the  bandwidth  versus  array  size  will  be  given. 
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Among  various  efforts  to  overcome  the  major  limitation  of  the  finite-element  time-domain  (FETD) 
method,  which  is  the  need  to  solve  a  global  matrix  equation  for  each  time  step,  one  important  progress  is 
the  development  of  the  dual-field  domain-decomposition  (DFDD)  method  (Z.  Lou  and  J.  M.  Jin,  “A  novel 
dual-field  time-domain  finite-element  domain-decomposition  method  for  computational 
electromagnetics,”  IEEE  Trans.  Antennas  Propag .,  vol.  54,  pp.  1850-1862,  June  2006).  In  DFDD,  the 
original  problem  is  broken  down  into  smaller  subdomain  problems  so  that  only  matrix  equations  at  the 
subdomain  level  are  solved.  In  contrast  to  FETD  which  solves  one  2nd-order  wave  equation,  DFDD  solves 
two  2nd-order  wave  equations: 


I 


— (VxT)-(VxE) 
Mr 


l 


!(VxT).(VxH)  +  4t.^ 

Sr  C()  dt 


dV  =  s,  j]  (nxT)|nx 


(1) 

(2) 


where  and  Ms  are  the  equivalent  electric  and  magnetic  currents  defined  at  the  subdomain  interfaces. 
Eqs.  (1)  and  (2)  are  solved  at  integer  and  half-integer  time  steps,  respectively,  and  hence  the  system  is 
marched  in  a  leapfrog  manner.  By  computing  the  equivalent  currents  as  J s  =nx Hneighbor and 

Ms  =  -n  x  Eneighbor  where  h  denotes  the  vector  normal  to  the  interface  and  pointing  outward  and  Eneigbbor 
and  Hneighbor  are  the  fields  from  the  neighboring  subdomain,  information  is  exchanged  among  subdomains 
and  the  tangential  field  continuities  are  weakly  enforced. 


The  DFDD  method  has  a  clear  physical  meaning.  Eq.  (1)  represents  the  application  of  the  surface 
equivalent  principle  with  a  surface  electric  current  radiating  in  front  of  a  magnetic  wall,  while  Eq.  (2) 
represents  the  same  principle  with  a  surface  magnetic  current  radiating  in  front  of  an  electric  wall.  Eqs. 
(1)  and  (2)  also  reflect  the  well-known  Huygens’  principle,  since  the  wave  front  at  the  subdomain 
interfaces  can  be  viewed  as  a  secondary  source  that  produces  the  field  inside  the  local  subdomain. 


Another  important  domain-decomposition  method  is  the  discontinuous  Galerkin  time -domain  (DGTD) 
method,  which  can  be  categorized  into  two  versions,  DGTD-Central  and  DGTD-Upwind,  depending  on 
whether  the  central  or  upwind  numerical  fluxes  are  used.  In  terms  of  the  interpolation  error,  DFDD, 
DGTD-Central,  and  DGTD-Upwind  have  the  convergence  rates  of  0(hp) ,  0(hp) ,  and  0(hp+l), 
respectively,  where  h  is  the  mesh  size  and  p  is  the  polynomial  order  of  the  basis  functions.  In  terms  of  the 
energy-conservation  property,  DFDD  and  DGTD-Upwind  have  no  numerical  dissipation  while  DGTD- 
Upwind  has  slight  numerical  dissipation  which  can  be  reduced  by  employing  higher-order  basis 
functions. 


With  respect  to  efficiency,  our  numerical  tests  show  that  the  CPU  time  for  DFDD  and  DGTD-Central  are 
about  the  same,  while  DGTD-Upwind  consumes  more  time  than  the  other  two  methods,  when  the  same 
mesh  and  maximum  allowable  time  step  are  used  for  all  three  methods.  When  the  mesh  is  fixed  and  the 
order  of  basis  functions  is  increased,  the  three  methods  consumes  about  the  same  amount  of  time  for  a 
specified  error  level.  When  the  order  of  basis  functions  is  fixed  and  the  mesh  is  refined,  DGTD-Upwind 
outperforms  the  other  two  methods  since  it  has  a  higher  error  convergence  rate. 
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In  this  talk,  we  will  present  a  new  full  vectorial  generalized  discontinuous  Galerkin  beam  propagation  method 
(GDG-  BPM)  to  accurately  handle  the  discontinuities  in  electromagnetic  fields  associated  with  wave  propagations 
in  inhomogeneous  optical  waveguides.  The  numerical  method  is  a  combination  of  the  traditional  beam  propagation 
method  (BPM)  with  a  newly  developed  generalized  discontinuous  Galerkin  (GDG)  method  .  The  GDG  method  is 
based  on  a  reformulation,  using  distributional  variables  to  account  for  solution  jumps  across 
material  interfaces,  of  Schro  'dinger  equations  resulting  from  paraxial  approximations  of  vector  Helmholtz  equations. 
Four  versions  of  the  GDG-BPM  are  obtained  for  either  the  electric  or  magnetic  field  components.  Modeling  of 
wave  propagations  in  various  optical  fibers  using  the  full  vectorial  GDG-BPM  is  included.  Numerical  results 
validate  the  high  order  accuracy  and  the  flexibility  of  the  method  for  various  types  of  interface  jump  conditions. 
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Complex  electromagnetic  applications  require  sophisticated  methods  to  solve  the  time-dependent 
Maxwell’s  equations.  Discontinuous  Galerkin  (DG)  finite  element  methods  are  highly  suitable 
for  time-domain  simulations.  They  support  various  types  and  shapes  of  elements,  non-conformal 
meshes  and  non-uniform  degrees  of  approximation.  Furthermore,  they  can  lead  into  fully  explicit 
time-marching  schemes.  Finally,  information  exchange  is  required  only  between  neighboring  el¬ 
ements,  which  leads  to  high  parallelization  efficiency.  In  electromagnetism,  DG  methods  were 
recently  applied  for  the  solution  of  the  time  dependent  Maxwell’s  equations.  In  this  paper  an  In¬ 
terior  Penalty  Discontinuous  Galerkin  (IPDG)  method  will  be  presented  to  solve  the  time  domain 
Maxwell’s  equations.  More  specifically,  emphasis  will  be  given  to  the  case  of  non-conformal  meshes. 
The  ability  to  handle  non-conformal  meshes  can  be  of  great  practical  importance.  In  many  cases 
the  high  level  of  complexity  found  in  practical  applications  (i.e.  IC  packaging  analysis)  makes 
the  task  of  generating  a  conformal  mesh  for  the  whole  structure  extremely  difficult.  Therefore,  a 
divide  and  conquer  approach  is  needed  in  such  cases.  In  our  approach,  the  computational  domain 
is  initially  decomposed  into  non-overlapping  sub-domains.  Currently,  only  planar  interfaces  are 
supported  between  subdomains.  Then,  each  sub-domain  is  meshed  completely  independently  re¬ 
sulting  in  non-conformal  meshes  at  subdomain  interfaces.  This  approach  provides  great  flexibility 
in  the  meshing  process.  Moreover,  to  take  advantage  of  the  high  parallelization  efficiency  offered  by 
discontinuous  Galerkin  methods  an  MPI  parallelization  has  been  implemented.  Each  sub-domain 
Mi,  is  mapped  to  an  MPI  process  and  each  MPI  process  is  associated  with  a  cluster  node.  Within 
each  node,  if  needed,  all  CPU  cores  are  used  with  OpenMP  to  update  the  sub-domain.  This  can 
be  considered  as  a  hybrid  MPI+OpenMP  approach.  Furthermore,  due  to  the  multi-scale  nature 
of  most  application  the  generated  meshes  often  contain  very  small  or  distorted  elements.  There¬ 
fore,  the  stability  condition  will  result  in  a  very  small  time  step  St.  To  guarantee  stability  for  all 
the  elements  we  must  choose  St  —  Stmin  —  min(Sti ),  i.e.  Stmin  is  the  minimum  of  all  the  local 
min(Sti).  Consequently,  the  CPU  time  will  significantly  increase.  Thus,  to  mitigate  this  problem, 
we  combine  our  parallelization  with  a  local  time-stepping  strategy  to  increase  efficiency  and  reduce 
the  computational  time. 
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Ideal  magnetohydrodynamics  (MHD)  equations  arise  in  many  areas  such  as  astrophysics,  space 
physics  and  engineering.  They  consist  of  a  set  of  nonlinear  hyperbolic  conservation  laws,  with  a 
divergence- free  constraint  on  the  magnetic  filed.  Neglecting  this  constraint  in  numerical  simulation 
may  lead  to  instability  or  non-physical  features  in  solutions.  Another  numerical  issue  which  can  lead 
to  numerical  instability  is  the  negative  pressue.  In  this  talk,  I  will  present  our  work  in  developing 
numerical  divergence  treatments  and  positivity  preserving  techniques  in  the  discontinuous  Galerkin 
(DG)  setting. 

In  (F.  Li  and  C.-W.  Shu,  Journal  of  Scientific  Computing,  v22-23  (2005),  pp. 413-442),  we  introduced 
a  locally  divergence-free  DG  method  for  the  ideal  MHD  equations.  Compared  with  the  standard 
DG  methods,  the  new  component  is  the  use  of  the  locally  divergence-free  approximations  for  the 
magnetic  field.  Though  the  idea  is  simple,  the  resulting  method  demonstrates  comparable  or 
better  accuracy  as  well  as  good  stability  in  numerical  simulations,  with  reduced  computational 
complexity.  On  the  other  hand,  the  divergence  of  the  numerical  magnetic  field  is  still  nonzero  due 
to  the  jump  in  its  normal  component  across  element  interfaces.  Motivated  by  this,  we  recently 
proposed  a  family  of  exactly  divergence-free  central  DG  methods  for  the  ideal  MHD  system  in  (F. 
Li,  L.  Xu  and  S.  Yakovlev,  Journal  of  Computational  Physics,  v230  (2011),  pp. 4828-4847)  and 
(F.  Li  and  L.  Xu,  Journal  of  Computational  Physics,  accepted  (2011)).  While  other  conservative 
quantities  are  evolved  with  standard  central  DG  methods,  the  magnetic  field  is  updated  such 
that  its  normal  component  is  first  approximated  along  the  mesh  interfaces  based  on  the  magnetic 
induction  equations.  A  divergence- free  reconstruction  is  then  carried  out  element  by  element.  To 
achieve  higher  order  accuracy  and  to  ensure  the  unique  solvability  of  the  reconstruction,  additional 
information  for  the  magnetic  field  is  extracted  by  revisiting  the  magnetic  induction  equation.  The 
overall  algorithm  is  local,  and  it  can  be  of  arbitrary  order  of  accuracy.  Numerical  examples  are 
presented  to  demonstrate  the  performance  of  these  methods. 

By  following  the  work  of  (Zhang  and  Shu,  Journal  of  Computational  Physics,  v229  (2010),  pp.8918- 
8934),  we  further  develop  the  positivity  preserving  DG  and  central  DG  methods  in  MHD  simu¬ 
lations.  This  is  based  on  a  simple  limiter  which  can  be  proved  to  maintain  uniform  high  order 
accuracy.  This  limiter  can  be  easily  implemented  in  multi-dimensions  and  on  general  meshes. 
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A  new  kind  of  prism  element  with  a  triangular  base  is  presented  for  discretization  of  layered  struc¬ 
tures.  For  the  triangular  base,  Mixed-order  Curl-Conforming  Edge-Based  Vector  basis  functions 
(CtLn  and  LtQn)  are  used  (A.  F.  Peterson,  S.  L.  Ray,  and  R.  Mittra,  Computational  methods  for 
electromagnetics ,  IEEE  press  New  York,  1998);  the  prism  height  adopts  spectral  basis  functions 
based  on  Gauss-Lobatto-Legendre  polynomials,  with  an  arbitrary  order  of  interpolation  (J.  Lee, 
T.  Xiao,  and  Q.  H.  Liu,  “A  3-D  Spectral-Element  Method  Using  Mixed-Order  Curl-Conforming 
Vector  Basis  Functions  for  Electromagnetic  Fields”,  IEEE  Transactions  On  Microwave  Theory  And 
Techniques,  Vol.  54,  No.  1,  Jan  2006).  These  two  families  of  basis  functions  allow  three  major  ad¬ 
vantages:  triangle  shapes  for  a  flexible  description  of  the  geometry  in  each  layer  plane;  a  high-order 
accuracy  and  a  low  memory  cost  in  the  height  direction;  and,  an  inherent  orthogonality  between 
basis  functions  in  the  layer  plane  (triangle)  and  height  (spectral). 

The  eigenvalues  and  eigenvectors  are  evaluated  with  this  element,  and  the  results  are  compared 
with  the  Finite  Elements  Method  (FEM)  based  on  tetrahedral  elements  and  Spectral  Elements 
Method  (SEM)  based  on  hexahedral  elements.  Excellent  accuracy  and  performance  of  the  new 
prism  element  are  shown  for  wave  equation  and  Maxwell’s  equations.  From  this  analysis,  it  is 
verified  that,  in  order  to  eliminate  spurious  solutions,  one  must  employ  mixed-order  of  interpolation 
for  electric  and  magnetic  fields  (J.  Chen,  Q.  H.  Liu,  M.  Chai,  and  J.  A.  Mix,  “A  Nonspurious  3- 
D  Vector  Discontinuous  Galerkin  Finite-Element  Time-Domain  Method”,  IEEE  Microwave  And 
Wireless  Components  Letters,  Vol.  20,  No.  1,  January  2010).  Results  in  time  domain  for  a  single 
domain  also  show  good  performance. 

An  obvious  application  of  this  new  element  is  layered  structures,  with  arbitrary  geometries  and 
scales  for  different  layers,  for  instance  on-chip  interconnects.  For  this  kind  of  applications,  in  this 
work  each  layer  is  considered  as  a  subdomain  with  independent  mesh,  orders  of  interpolations,  and 
shared  boundaries  with  neighbor  subdomains.  To  solve  the  connection  between  sudbomains  the 
discontinuous  Galerkin  method  with  a  Riemann  solver  is  applied.  Furthermore,  subdomains  are 
assumed  placed  in  a  sequential  order  (very  common  for  layered  structures)  and  consequently  the 
coupled  system  will  take  the  form  of  a  block  tri-diagonal  system;  thus,  a  block  Thomas  algorithm  is 
integrated  to  accelerate  and  improve  the  accuracy  of  the  implicit  Crank-Nicolson  method  (J.  Chen, 
L.  E.  Tobon,  M.  Chai,  J.  A.  Mix,  and  Q.  H.  Liu  “Efficient  Implicit-Explicit  Time  Stepping  Scheme 
With  Domain  Decomposition  for  Multiscale  Modeling  of  Layered  Structures” ,  IEEE  Transactions 
On  Components,  Packaging  And  Manufacturing  Technology,  Vol.  1,  No.  9,  September  2011). 
Finally,  some  applications  for  integrate  circuits  are  presented,  with  comparisons  with  traditional 
methods  and  commercial  software. 
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The  micro-resistivity  imager  is  a  very  useful  tool  for  oil  drilling  and  exploration.  During  rotating  /  sliding 
in  borehole,  it  can  generate  an  image  containing  facture  information  on  the  borehole  wall.  To  model  the 
micro-resistivity  imager  is  usually  a  multiscale  problem:  the  length  of  an  imager  tool  is  several  meters, 
while  the  details  on  the  imager  sensors  and  the  factures  on  the  borehole  wall  can  be  as  small  as 
millimeters  or  even  micrometers. 

In  our  study  we  need  time  domain  multiscale  simulation  of  micro-resistivity  imagers.  This  task  is  very 
challenging  for  conventional  methods,  such  as  the  popular  FDTD  (finite-element  time-domain)  technique. 
On  one  hand,  spatial  discretization  of  a  multiscale  structure  by  conventional  methods  oftentimes  leads  to  a 
huge  number  of  unknowns;  on  the  other  hand,  because  of  fine  cells  existing  in  the  discretized  multiscale 
system,  the  stability  criterion  in  time  stepping  can  be  extremely  strict,  and  an  unaffordable  number  of  time 
steps  is  required  to  finish  the  simulation  of  a  given  time  window. 

Here  we  employ  an  efficient  DG-FETD  (discontinuous  Galerkin  finite-element  time-domain)  method  for 
solving  this  multiscale  problem.  An  imager  structure  is  divided  into  several  subdomains  based  on  the 
electrical  size  of  each  part.  Large  and  higher  order  finite  elements  are  used  in  discretizing  electrically 
coarse  structures;  small  and  lower  order  finite  elements  are  used  to  capture  the  geometric  characteristics 
of  fine  structures.  The  Riemann  solver  is  employed  to  stitch  the  separate  subdomains  together.  After 
spatial  discretization,  the  implicit-explicit  Runge-Kutta  (IMEX-RK)  scheme  is  utilized  to  perform  time 
integration.  The  explicit  RK  scheme  is  used  for  electrically  coarse  subdomains  discretized  with  coarse 
mesh,  and  the  implicit  RK  scheme  is  used  to  overcome  the  stability  criterions  of  electrically  fine 
subdomains.  Numerical  examples  demonstrate  that  the  DG-FETD  method  can  be  several  orders  more 
efficient  than  conventional  methods  in  modeling  the  micro-resistivity  imager. 
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The  Improved  Concurrent  Electromagnetic  Particle-In-Cell  (ICEPIC)  code  is  used  to  model  high  power 
microwave  (HPM)  devices.  ICEPIC  is  developed  since  1994,  and  the  code  is  massively  parallel  with 
scalability  tested  to  over  16000  CPUs.  The  relevant  physics  to  model  HPM  devices,  including  particle 
emission,  relativistic  particle  tracking,  perfectly  matched  layer  absorbing  boundary  conditions,  and 
material  handling  are  included  in  the  code.  Temporal  filtering  is  used  to  remove  poorly  resolved  waves 
that  can  lead  to  non-physical,  numerical  Cerenkov  radiation.  The  code  also  contains  an  interface  to 
allow  other  FEM  codes  to  model  antenna  structures  with  additional  flexibility.  For  high  density  plasmas, 
the  PIC  method  becomes  intractable,  and  a  simplified  fluid  model  is  included  in  the  code.  The  model 
handles  both  unmagnetized  and  magnetized  plasmas. 

Recently,  there  is  interest  in  using  frequency  dispersive  and  non-linear  dielectric  materials  in  HPM 
devices.  The  code  includes  a  Lorentz-Drude  model  in  the  time  domain  for  inclusion  of  frequency 
dispersive  dielectrics.  Nonlinear  material  modeling  is  an  area  of  active  research.  Two  types  of  materials 
are  considered.  The  first  is  nonlinear  dielectric  materials  for  which  the  permittivity  is  a  function  of  the 
electric  field  magnitude.  The  second  is  ferrite  material.  For  these  materials,  the  Landau-Lifshitz-Gilbert 
equation  is  discretized.  Nonlinear  materials  can  suffer  from  a  non-physical  instability  that  is  related  to  the 
spatial  cutoff  frequency  of  the  computational  mesh.  Efforts  to  suppress  the  instability  by  filtering  are 
discussed. 
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The  most  straightforward  method  of  solving  nonlinear  electromagnetic  field  problems  in  the  time  domain 
by  the  method  of  finite  elements  is  using  time-stepping  techniques.  This  requires  the  solution  of  a  large 
nonlinear  equation  system  at  each  time  step  and  is  therefore  very  time  consuming,  especially  if  a  three- 
dimensional  problem  is  being  treated.  If  the  excitations  are  non-periodic  or  if,  in  case  of  periodic 
excitations,  the  transient  solution  is  required,  one  cannot  avoid  time-stepping.  In  many  cases  however,  the 
excitations  of  the  problem  are  periodic,  and  it  is  only  the  steady-state  periodic  solution  which  is  needed. 
Then,  it  is  wasteful  to  step  through  several  periods  to  achieve  this  by  the  “brute  force”  method  of  time 
stepping. 

A  time  domain  technique  using  the  fixed-point  method  to  decouple  the  time  steps  has  been  introduced  in 
(O.  Biro,  K.  Preis,  An  efficient  time  domain  method  for  nonlinear  periodic  eddy  current  problems,  IEEE 
Transactions  on  Magnetics,  42,  (2006),  695-698).  In  later  work,  the  optimal  choice  of  the  fixed  point 
reluctivity  for  such  problems  has  been  determined  both  in  the  time  domain  and  using  harmonic  balance 
principles.  The  method  can  be  applied  to  three-dimensional  problems  both  in  terms  of  a  magnetic  vector 
potential  and  an  electric  scalar  potential  (A,v-A  formulation)  and  applying  a  current  vector  potential  and  a 
magnetic  scalar  potential  (T,0-0  formulation). 

The  aim  of  this  paper  is  to  show  the  application  of  the  method  to  industrial  problems  arising  in  the  design 
of  large  power  transformers. 
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Hybrid  FETD/FDTD  techniques  that  combine  the  advantages  of  both  finite-element 
time-domain  (FETD)  and  finite-difference  time-domain  (FDTD)  methods  and  avoid 
their  weaknesses  have  been  developed  recently.  The  FETD  method  is  applied  to 
model  structures  with  complex  geometries  and  boundaries,  and  the  FDTD  method  is 
efficient  to  model  large  homogeneous  volumes  and  more  regular  structures.  In  the 
previous  hybrid  approaches,  to  our  knowledge,  the  structured  FDTD  grid  and 
unstructured  FETD  mesh  are  required  to  be  conformal  to  each  other;  thus,  a  transition 
element  is  needed  to  joint  the  structured  FDTD  grid  with  the  unstructured  FETD  mesh, 
making  mesh  generation  very  challenging.  To  improve  flexibility  and  efficiency  of  the 
FETD/FDTD  hybridization,  a  non-conforming  mesh  between  the  structured  and 
unstructured  grids  is  desirable  for  simulating  multi-scale  structures.  We  recently 
proposed  the  highly  parallel  discontinuous  Galerkin  method  as  an  appropriate 
framework  for  hybridizing  the  FETD  and  FDTD  methods,  as  the  discontinuous 
Galerkin  method  allows  discontinuous  basis  functions  across  different  subdomains. 

In  this  presentation,  we  summarize  the  hybrid  FETD/FDTD  technique  based  on  the 
discontinuous  Galerkin  method  developed  recently.  This  technique  can  use  arbitrary 
non-conformal  meshes  between  the  structured  FETD  mesh  and  unstructured  FDTD 
grid.  Multiple  FETD  subdomains  are  allowed  in  this  method,  so  it  is  suitable  to  be 
applied  to  multiple  zones  divided  for  parallel  computation.  In  the  hybrid  method, 
either  explicit  scheme  such  as  the  4th-order  Runge-Kutta  method  (RK4)  or  implicit 
scheme  such  as  Crank-Nicholson  method  can  be  used  in  the  FETD  regions  for  time 
integration.  Several  examples  will  be  presented  to  show  the  validity  and  efficiency  of 
the  proposed  hybrid  method. 
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In  this  paper,  we  present  our  efforts  in  combating  the  multi-scale  thermal  problems  through  the  use  of 
non-conformal  Domain  Decomposition  Method  (DDM).  The  non-conformal  DDM  have  received 
considerable  attention  because  they  provide  effective,  efficient  preconditioned  iterative  solution 
algorithms  for  the  numerical  solution  of  moderately  stiff  parabolic  partial  differential  equations  (PDEs)  in 
inhomogeneous  domains. 

The  numerical  procedure  that  we  introduced  herein  is  based  on  the  non-conformal  and  non-overlapping 
DDM  for  transient  thermal  analysis.  Specifically,  we  partition  the  problem  geometry  into  many  non¬ 
overlapping  sub-domains.  Within  each  sub-domain,  we  strive  to  maintain  similar  geometrical  feature 
sizes.  Next,  each  of  the  sub-domains  is  meshed  independently  without  considerations  of  adjacent  sub- 
domains.  The  required  temperature  and  heat  flux  continuities  across  domain  interfaces  are  guaranteed  by 
the  Robin  transmission  condition,  which  is  enforced  weakly  through  the  Galerkin  method.  Consequently, 
the  proposed  non-conformal  DDM  greatly  relaxes  the  burden  of  mesh  generation  for  complicated  multi¬ 
scale  thermal  problems  and  provides  an  inherent  parallelism  and  flexibility. 

Moreover,  the  time  discrelization  employed  is  based  on  an  unconditional  stable,  implicit  Euler  scheme. 
The  unconditionally  stable  time  marching  algorithm  is  beneficial  since  the  time-step  size  is  no  longer 
governed  by  the  spatial  discretization  of  the  mesh,  but  rather  by  the  desired  accuracy.  Additionally,  the 
numerical  investigations  of  the  convergence  properties  of  the  proposed  non-conformal  DDM  are 
included. 


Furthermore,  we  extend  the  non-conformal  DDM  formulation  for  the  steady  state  thermal  analysis  to 
study  the  steady  thermal  temperature  distribution  within  complex  ICs.  Finally,  numerical  results  are 
shown  for  a  chip-package-PCB  example  with  cooling  techniques  of  both  natural  convection  and  forced 
convection  cooling  through  the  heat  sinks. 
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Time-domain  methods  can  be  categorized  into  two  broad  classes.  One  is  the  explicit  time-domain 
method;  the  other  is  the  implicit  time-domain  method.  Explicit  methods  can  avoid  solving  a  matrix,  while 
implicit  methods  require  a  matrix  solution.  Despite  its  advantage  of  being  matrix  free,  an  explicit  method 
requires  the  time  step  to  be  restricted  by  the  smallest  space  step  for  ensuring  stability.  For  problems  that 
have  fine  features  relative  to  working  wavelength  like  integrated  circuits,  explicit  methods  require  a  large 
number  of  time  steps  to  finish  one  simulation,  which  is  computationally  expensive.  Existing 
unconditionally  stable  methods  (methods  that  permit  the  use  of  any  large  time  step  without  becoming 
unstable)  are  all  implicit  methods  such  as  ADI  (alternating  direction  implicit)-FDTD,  the  CN  (Crank- 
Nicolson)-FDTD,  the  Laguerre-FDTD,  and  the  Newmark-Beta  based  time -domain  finite-element  method 
(TDFEM).  Although  they  permit  the  use  of  any  large  time  step,  they  require  the  solution  of  a  large-scale 
matrix  for  large  problems.  Moreover,  late-time  instability  has  been  observed  from  implicit  methods. 
Recently,  research  has  also  been  done  to  extend  the  stability  limit  of  the  explicit  FDTD  method  by  spatial 
filtering.  As  yet,  no  explicit  methods  have  achieved  unconditional  stability.  The  research  question 
considered  in  this  work  is:  can  an  explicit  method  be  made  unconditionally  stable  so  that  its  matrix-free 
strength  can  be  retained  and  its  shortcoming  in  time  step  can  be  eliminated? 

The  preliminary  research  of  this  work  was  reported  in  (Q.  He  and  D.  Jiao,  2011  IEEE  AP-S 
symposium),  where  we  quantitatively  identified  the  root  cause  of  the  instability  associated  with  an 
explicit  time-domain  method  that  uses  a  time  step  beyond  stability  criterion.  Based  on  such  a  finding,  we 
successfully  developed  an  explicit  time -domain  finite-element  method  that  is  unconditionally  stable,  i.e., 
stable  for  any  large  time  step.  The  formulations  given  in  (Q.  He  and  D.  Jiao,  2011  IEEE  AP-S 
symposium)  are  for  lossless  cases  in  which  only  ideal  dielectrics  and  conductors  are  considered.  In 
addition,  linear  (optimal)  computational  complexity  was  not  reported.  The  contribution  of  this  paper  is  the 
successful  development  of  a  linear-complexity  explicit  and  unconditionally  stable  time-domain  finite- 
element  method  for  the  simulation  of  arbitrarily  shaped  3-D  non-ideal  conductors  embedded  in 
inhomogeneous  materials.  We  also  provide  a  theoretical  basis  for  making  an  explicit  time-domain  method 
unconditionally  stable.  Numerical  experiments  have  demonstrated  that  the  proposed  unconditionally 
stable  explicit  method  outperforms  both  the  conditionally  stable  explicit  method  and  the  unconditionally 
stable  implicit  method  in  computational  efficiency  and  stability.  Although  the  proposed  method  is 
presented  in  the  framework  of  a  time-domain  finite-element  method,  the  essential  idea  can  be  applied  to 
other  time-domain  methods. 
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Discontinuous  Galerkin  (DG)  finite  element  method  is  well-suited  on  unstructured  meshes  for  high 
order  approximation  with  the  freedom  of  non-matching  grids  and  choosing  the  order  of  basis  func¬ 
tions  in  each  element  locally.  Moreover,  DG  can  handle  complicated  geometries  with  curved  bound¬ 
ary  easily.  Further,  since  information  exchange  in  DG  only  involves  neighboring  elements  which 
results  in  block  diagonal  mass  matrices,  high  efficiency  in  parallelization  can  be  achieved  for  time 
domain  method  using  explicit  time  marching  scheme.  As  a  truncation  to  the  computational  do¬ 
main  for  finite  element  method, perfect  matched  layer  (PML)  is  adopted  in  this  paper.  By  making 
the  shape  of  the  layer  conformal  to  the  target  scatterer  the  size  of  the  domain  can  be  reduced. 
Conformal  PML  can  be  treated  as  an  artificial  anisotropic  layer  where  the  material  properties  vary 
continuously  in  space.  In  order  to  model  this  spatially  changing  materials  in  conformal  PML,  this 
paper  uses  a  universal  array  approach  where  inside  each  element  the  material  property  tensors 
change  as  a  function  of  location.  Compared  to  the  conventionally  treatment  of  element-wise  con¬ 
stant  material  property,  this  method  provides  a  more  accurate  physical  model  for  problems  with 
continuously  varying  material  properties  in  an  efficient  way,  in  this  case  the  modeling  of  conformal 
PML  and  many  other  real  life  applications  as  well. 

The  main  contribution  of  this  paper  is  a  combination  of  the  universal  array  method  for  continuous 
material  property  and  MPI-based  parallel  implementation  of  three  dimensional  (3D)  DGTD  method 
with  explicit  time  marching  scheme  and  local  time  stepping  to  solve  for  the  first  order  Maxwell’s 
equations. 
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Hybrid  finite  element  method  and  method  of  moments  is  a  powerful  tool  to  simulate  electromagnetic 
scattering  of  dielectric  object.  The  finite  element  method  can  model  complicated  structure  with 
unstructured  mesh.  The  method  of  moments  exactly  describe  the  propagation  of  electromagnetic 
field  in  the  background.  However,  when  the  size  of  the  simulated  system  is  much  larger  than  the 
wavelength,  it  is  more  difficult  for  this  solver  to  find  the  solution.  When  the  number  of  unknowns 
become  larger,  the  solver  consume  more  memory.  In  addition,  the  iterative  solver  might  not  be 
able  to  converge  when  the  size  of  the  system  matrix  is  large. 

One  of  the  solution  to  this  bottleneck  is  to  introduce  domain  decomposition,  which  split  up  the 
whole  system  into  a  few  smaller  system,  or  sub-domains.  The  system  matrix  of  each  sub-domain 
is  smaller  enough  to  be  solved  by  directed  solver.  The  solution  of  the  whole  system  is  obtained 
by  considering  the  interfere  between  the  solution  of  adjacent  sub-domains.  We  use  Riemann  solver 
to  model  the  interfere  between  adjacent  sub-domains.  The  electromagnetic  field  within  each  sub- 
domain  is  model  to  spectral  element  method,  which  is  a  kind  of  high  order  finite  element  method 
base  on  hexahedral  elements.  The  Riemann  solver  is  inserted  into  the  surface  integral  of  the  weak 
form.  For  finite  periodic  system,  each  period  is  defined  as  a  sub-domain.  Thus,  each  sub-domain 
have  the  same  structure,  which  generates  the  same  system  matrix.  As  a  result,  we  don’t  need  to 
store  the  repeating  system  matrix. 

The  background  is  defined  as  one  sub-domain,  which  is  modeled  by  method  of  moments.  In  order 
to  reduce  memory  cost  and  CPU  time,  the  spectral  integral  method  is  used  to  store  and  solve  the 
method  of  moments  matrixes.  The  combined  field  integral  equation  is  used  in  the  implementation 
of  the  method  of  moments.  Numerical  examples  show  that  the  convergent  speed  is  fast. 
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The  use  of  vector  finite  elements  for  representing  electromagnetic  fields  is  standard  procedure 
for  modern  simulations  and  is  a  staple  of  academic  instruction.  Properly  constructed  edge 
elements  offer  clear  mathematical  advantages  such  as  tangential  field  component  continuity,  the 
suppression  of  spurious  solutions,  and  the  satisfaction  of  the  Nedelec  conditions  to  prevent 
pathological  representations.  Two  penalties  that  are  tolerated  to  obtain  the  advantages  of  edge 
elements  are  their  mixed  polynomial  order  (edge  vs.  normal  field  descriptions)  and  the 
infrastructural  overhead  associated  with  a  non-nodal  mesh  bookkeeping.  In  the  first  case,  the 
same  mathematical  advantages  that  allow  vector  bases  to  simultaneously  satisfy  tangential 
boundary  conditions,  maintain  consistency  with  the  operator,  and  retain  the  physics  of  the 
problem  result  in  the  tangential  and  normal  field  components  at  the  edges  to  be  of  mixed 
polynomial  order.  While  many  elegant  forms  have  been  published  to  make  edge  elements 
hierarchical  to  increased  polynomial  order,  none  overcome  the  obstacle  that  the  edge  field 
description  is  of  order  n  -  1  compared  to  polynomial  order  n  interior  to  the  triangle.  The  second 
issue  of  not  using  a  straightforward  nodal  representation  common  to  scalar  field  problems  is  a 
practical  one.  For  example,  even  order  edge  field  basis  functions  (such  as  published  constant  or 
quadratic  tangential  field  elements)  require  not  only  spatial  coordination  between  adjacent 
triangles,  but  edge  directionality  to  be  consistent.  This  necessitates  cumbersome  modifications  or 
add-ons  to  the  prevalent  technology  infrastructure  for  mesh  generation.  Custom  software 
requirements  are  a  serious  hindrances  to  advancing  the  using  of  more  sophisticated  basis 
functions  in  electromagnetic  simulations. 

In  this  work,  we  propose  alternative  polynomial  basis  functions  for  use  in  electromagnetic  field 
calculations  that  can  be  used  with  triangular  elements  on  a  simple  nodal  mesh.  We  employ 
scalar,  fifth-order  Hermite  interpolation  polynomials  to  solve  Maxwell’s  equations  in  two 
dimensions  in  the  finite  element  method.  Each  field  component  is  represented  by  scalar  Hermite 
elements  while  maintaining  field  and  derivative  continuity  for  treating  boundary  conditions  and 
avoiding  spurious  modes.  We  have  analyzed  homogeneous  conducting  waveguides, 
inhomogeneous  interior  waveguides,  and  photonic  crystals.  The  Hermite  interpolation  functions 
provide  greater  accuracy  than  vector  finite  elements  of  equal  polynomial  order,  while  bypassing 
the  issue  of  spurious  modes  observed  when  using  Lagrange  polynomials.  Comparisons  are  made 
to  edge  vector  element  calculations  with  comparable  number  of  degrees  of  freedom.  In  addition 
to  highly  accurate  eigenvalue  computation,  the  calculated  spatial  fields  representing  the 
eigenmodes  do  not  suffer  degradation  from  mixed  polynomial  order.  The  impact  is  that  few 
elements  can  be  used  to  resolve  the  same  level  of  accuracy. 

*Approved  for  Public  Release,  Distribution  Unlimited 
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Parallel  Finite  Element  Electromagnetic  Code  Suite  for  Accelerator  Modeling 

and  Simulation 

Cho  Ng,  Lixin  Ge,  Kwok  Ko,  Kihwan  Lee,  Zenghai  Li,  Greg  Schussman,  and  Liling  Xiao 

SLAC  National  Accelerator  Laboratory 


A  parallel  scalable  electromagnetic  code  suite,  ACE3P,  has  been  developed  at  SLAC  over  the  past 
decades  for  particle  accelerator  design,  optimization  and  analysis.  Based  on  the  high-order  finite-element 
method,  ACE3P  allows  high-fidelity  modeling  of  complex  accelerator  structures  and,  running  on  state- 
of-the-art  supercomputers,  enables  computationally  challenging,  large-scale  problems  to  be  solved  with 
unprecedented  accuracy.  ACE3P’s  comprehensive  simulation  capabilities  include  eigensolvers  and 
scattering  matrix  calculation  in  the  frequency  domain,  transient  and  full-wave  particle-in-cell  simualtions 
in  the  time  domain,  charged  particle  tracking,  as  well  as  integrated  multiphysics  analysis  with 
electromagnetic,  thermal  and  mechanical  effects.  Applications  of  ACE3P  to  accelerator  projects  for 
wakefield  computation  at  the  system  level,  large-scale  optimization  of  accelerator  structures,  and  far-field 
radiation  calculation  for  photonic  bandgap  optical  fibers  will  be  presented. 
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Optical,  Electromagnetics,  and  Thermal  Modeling  of  Interaction  of  a  Focused 
Beam  of  Light  with  Plasmonic  Nanoparticles 

Eren  S.  Unlu1  and  Kursat  Sendur*1 

1  Sabanci  University,  Orhanli  -  Tuzla,  34956,  Istanbul,  Turkey,  sendur@sabanciuniv.edu 


The  interaction  of  photons  with  metallic  nanoparticles  is  important  to  a  number  of  emerging 
nanotechnology  applications  due  to  the  large  enhancement  and  tight  localization  of 
electromagnetic  fields  in  the  vicinity  of  nanoparticles.  This  interaction  has  existing  and  emerging 
applications  at  the  nanoscale,  including  near-field  scanning  optical  microscopy,  high-density  data 
storage,  and  bio-chemical  sensing.  Although  there  has  been  much  effort  to  understand  the  effects 
of  various  parameters  on  the  plasmon  resonances  of  nanoparticles,  the  interaction  of 
nanoparticles  with  a  focused  beam  of  light  has  largely  been  omitted  in  the  context  of  particle 
plasmons.  In  addition,  a  tightly  focused  beam  of  incident  light  can  be  used  to  excite  a  metallic 
nano-antenna.  Since  the  incident  electric  field  illuminating  the  nano-antenna  increases,  the  near¬ 
field  radiation  from  the  nano-antenna  also  increases. 


In  the  first  part  of  this  study,  it  is  demonstrated  that  the  focused  light  can  be  utilized  to 
manipulate  the  distribution  of  Mie-plasmons  on  nanoparticles.  This  interaction  is  particularly 
important  for  applications  such  as  optical  and  magneto-optical  data  storage,  where  strong  optical 
fields  are  desired  at  specific  locations.  The  interaction  of  light  with  a  nanoparticle  is  investigated 
for  beams  with  various  angular  spectra.  To  analyze  the  effect  of  the  angular  spectrum  on  the 
particle  plasmons,  a  silver  nanoparticle  is  illuminated  using  a  focused  beam  of  light  with  various 
half-beam  angles.  Also  in  this  part,  interaction  of  a  focused  beam  of  light  with  a  nano-antenna 
structure  is  studied.  A  metallic  antenna  is  placed  in  the  vicinity  of  a  solid  immersion  lens,  and  the 
effect  of  various  parameters,  including  the  angular  spectrum  of  the  beam  on  the  near-field 
radiation  from  the  nano-antenna  is  studied. 


In  the  second  part  of  this  study,  localized  radiative  energy  transfer  from  a  near-field  emitter  to  a 
magnetic  thin  film  structure  is  investigated.  A  magnetic  thin  film  stack  is  placed  in  the  near-field 
of  the  plasmonic  nanoantenna  to  utilize  the  evanescent  mode  coupling  between  the  nanoantenna 
and  magnetic  thin  film  stack.  A  bow-tie  nano-optical  antenna  is  excited  with  a  tightly  focused 
beam  of  light  to  improve  near-field  radiative  energy  transfer  from  the  antenna  to  the  magnetic 
thin  film  structure.  A  tightly  focused  incident  optical  beam  with  a  wide  angular  spectrum  is 
formulated  using  Richards-Wolf  vector  field  equations.  Localized  radiative  energy  transfer 
between  the  near-field  emitter  and  the  magnetic  thin  film  structure  is  quantified  for  a  given 
optical  laser  power  at  various  distances  between  the  near-field  emitter  and  magnetic  thin  film. 
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Finite  element  analysis  of  wideband  nanostructures 
for  photovoltaic  applications 

Niccolo  Breda,  Giuseppe  Pelosi,  Stefano  Selleri,  and  Ruggero  Taddei 
Department  of  Electronics  and  Telecommunications,  University  of  Florence, 

Via  di  S.  Marta  3,  50139,  Florence,  Italy 

Nanostructures  enhanced  solar  cells  are  progressively  emerging  as  a  theoretically  viable  alternative  to  con¬ 
ventional  ones  (K.  R.  Catchpole  and  A.  Polman,  ’’Plasmonic  solar  cells”,  Optical  Express,  Vol.  16,  2008,  pp. 
21793-21800).  While  first  and  second  generations  solar  cells  exploit  traditional  photovoltaic  effect,  the  so 
called  third  generation  of  solar  energy  harvesting  techniques  is  aimed  at  taking  advantage  of  the  plasmonic 
resonance  phenomenon. 

Focalised  surface  plasmon  (ESP)  resonances  have  been  shown  to  obtain  a  great  electric  field  enhancement 
in  the  vicinity  of  metallic  nanoparticles  immersed  in  a  conventional  semi-conductive  slab  (Crozier,  K.B  et 
al.,  ’’Field  enhancement  and  resonance  in  optical  antennas”,  Quantum  Electronics  and  Laser  Science,  2003. 
QELS.  Postconference  Digest,  2003).  This  increase  in  the  electric  field  magnitude  can  then  be  exploited  to 
lead  to  greater  efficiencies  with  respect  to  conventional  solar  cells. 

The  LSP  resonance  occurs  at  a  particular  frequency  in  the  optical  range,  where  metals  cannot  be  modelled 
any  more  as  perfect  conductors,  but  rather  assume  a  plasma  behavior,  expressed  through  a  modified  Drude- 
Lorentz  model  (B.Ung,  Y.  Sheng,  ’’Interference  of  surface  waves  in  a  metallic  nanoslit”,  Optics  Express, 
Vol.  15,  No.  3,  2007,  pp.  1182-1190). 

To  gain  more  insight  on  this  particular  behavior,  finite  element  (FE)  simulations  are  conducted  on  a  periodic 
arrangement  of  gold  and  silver  nanoparticles  on  top  of  a  silicon  slab.  Particles  with  various  geometries 
are  investigated;  their  size  is  shown  to  affect  the  onset  of  the  plasmonic  resonance  phenomenon.  Particles 
with  sizes  between  20  and  50  nm  are  shown  to  lead  to  the  desired  field  enhancement,  which  appears  as  a 
minimum  in  the  reflection  coefficient.  The  analysis  is  then  deepened  by  placing  an  arrangement  of  metal 
particles  with  different  sizes  on  a  same  substrate. 
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Recently,  metal  and  carbon  nano  wires  (NWs)  have  received  a  serious  amount  of  interest  due  to 
their  potential  in  confining  light  transversally  to  sub-wavelength  dimensions  and  yet  to  be  used 
as  an  optical  waveguide  or  an  antenna  in  the  visible  (Bharadwaj  et  al .,  Adv.  Opt.  Photon.  1, 
2009).  Oulton  et  al.  has  shown  that  a  hybrid  optical  waveguide  (HOW),  which  consists  of  a 
dielectric  NW  separated  from  a  metal  surface  by  a  nano  scale  dielectric  gap,  can  provide  extremely 
long  propagation  length  (dozens  of  wavelengths)  and  strong  mode  confinement  (Nat.  Photonics 
2,  2008).  Experimental  and  theoretical  results  reveal  a  huge  potential  for  realistic  nano  scale 
semiconductor-based  plasmonics  and  photonics. 
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Figure  1:  Total  electric  field  (incident  +  scattered) 
for  2  NWs  (depicted  as  magenta  dashed  lines)  embed¬ 
ded  in  a  2-layered  medium.  Black  line  represents  the 
air-glass  interface.  White  lines  show  the  FEM  mesh. 
eNm  =  16>  eNW2  =  —7.6321  +  0.7306*. 


In  this  work  we  propose  a  novel  HOW 
structure  which  can  provide  even  longer 
propagation  lengths.  The  superiority 
of  the  design  comes  from  the  metallic 
loss  reduction  obtained  by  replacing  the 
metal  planar  layer  with  a  metal  cylin¬ 
der.  Numerical  results  (obtained  with  a 
two  dimensional  mode  solver)  show  that 
this  reduction  in  the  metal  surface  area 
increases  the  propagation  length. 

In  order  to  understand  electromagnetic 
field  distribution  throughout  the  whole 
structure,  we  also  calculated  the  total 
electric  field  due  to  plane  wave  inci¬ 
dence  of  different  incidence  angles.  In 
this  direction,  we  use  the  spectral  in¬ 
tegral  method  (SIM)  as  an  exact  ra¬ 
diation  boundary  condition  to  truncate 
the  computational  domain  in  the  finite 
element  method  (FEM)  to  form  a  hy¬ 
brid  SIM/FEM  which  is  applicable  to 
arbitrary  inhomogeneous  objects.  Fig. 
1  shows  the  field  distribution  around  2 
NWs,  where  the  silver  one  is  embedded 
in  glass  for  a  normally  incident  plane 
wave  (A  —  500  nm). 

At  the  conference,  we  will  discuss  the 
performances  of  the  proposed  and  tra¬ 
ditional  HOWs  in  great  detail. 
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Reduced-order  Electromagnetic  Modeling  in  the  Presence  of  Uncertainty 
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This  paper  is  concerned  with  the  development  of  reduced-order,  compact  models  that  account  for 
uncertainty  (randomness)  in  geometric  and  material  properties  in  the  finite  element  modeling  of 
electromagnetic  wave  phenomena  in  complex  structures.  Of  particular  interest  are  structures  where  such 
uncertainty  is  localized,  occurring  in  a  finite  number  of  regions  inside  the  structure.  For  such  cases,  the 
Monte  Carlo-based  numerical  electromagnetic  analysis  of  the  structure  can  be  expedited  by  eliminating 
the  repeated  mesh  generation  required  for  each  sampling  point  in  the  random  space  defined  by  the 
independent  random  variables  that  parameterize  the  uncertainty  in  the  structure. 

Toward  this  objective,  this  paper  considers  the  development  of  stochastic  reduced-order  macromodels  for 
each  one  of  the  regions  exhibiting  uncertainty,  with  each  macromodel  defined  over  a  fixed  surface 
enclosing  the  associated  random  region.  In  this  manner,  the  geometry  exterior  to  the  union  of  the  random 
regions  is  deterministically  defined.  Thus,  a  single  finite  element  mesh  needs  to  be  generated  for  the 
discretization  of  the  overall  structure.  The  only  thing  that  changes  for  each  sample  in  the  Monte  Carlo 
process  is  the  stochastic  macromodel  for  each  one  of  the  random  regions. 

The  stochastic  macromodel  is  defined  on  the  fixed  boundary  enclosing  the  random  region  in  terms  of  a 
global  admittance  matrix.  More  specifically,  with  the  tangential  electric  and  magnetic  fields  on  the 
boundary  interpolated  in  terms  of  appropriate  sets  of  deterministic  spatial  expansion  functions,  the  global 
admittance  matrix  relates  the  coefficients  in  these  expansions.  These  coefficients  are  functions  of  both 
frequency  and  the  independent  random  variables  involved  in  the  parameterization  of  the  uncertainty 
inside  the  region.  With  each  random  variable  defined  in  terms  of  its  probability  density  function,  use  is 
made  of  generalized  polynomial  chaos  expansions  for  the  representation  of  these  coefficients  in  terms  of 
truncated  polynomial  expansions  in  random  space.  With  a  finite  element  model  used  for  the  discretization 
of  the  interior  of  the  random  region,  a  sparse  stochastic  collocation  method  is  used  to  develop  reduced- 
order  macromodels  of  the  admittance  matrix  that  are  expressed  in  terms  of  generalized  polynomial  chaos 
expansions  over  the  random  space  for  the  region  (P.  Sumant,  et  al,  “A  Sparse  Grid  Based  Collocation 
Method  for  Model  Order  Reduction  of  Finite  Element  Approximations  of  Passive  Electromagnetic 
Devices  Under  Uncertainty,”  2010  IEEE  MTT-S  International  Microwave  Symposium  Digest,  pp.  1652  - 
1655,  Anaheim,  CA,  2010). 

The  basic  steps  involved  in  the  development  of  reduced-order,  stochastic  macromodel  is  presented  in  the 
context  of  specific  applications  involving  electromagnetic  wave  scattering  by  clusters  of  targets  exhibiting 
uncertainty.  These  applications  are  used  to  highlight  the  attributes  of  the  proposed  stochastic  finite 
element  macromodels  and  demonstrate  the  computational  benefits  resulting  from  their  utilization. 
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This  contribution  addresses  the  problem  of  applying  a  model  order  reduction  technique  locally  in 
three-dimensional  finite  element  analysis  of  waveguide  structures  and  resonators.  The  basic  idea 
of  the  presented  approach  is  to  identify  subregions  requiring  strong  mesh  refinement  and  replace 
their  large  subsystems  of  equations  with  the  much  smaller  ones,  called  macro-elements,  which  are 
created  by  means  of  the  model  order  reduction  ENOR  algorithm.  As  a  result,  the  overall  size  of 
the  problem  and  the  computation  time  are  significantly  reduced,  without  sacrificing  the  accuracy 
of  simulation.  What  is  more,  if  the  number  of  macro-elements  and  the  reduction  order  is  large, 
the  dense  macro-element  matrices  can  be  subject  to  the  second  reduction  and  the  diagonalization 
process. 

The  efficiency  of  the  reduction  process  depends  on  the  number  of  variables  inside  the  macro-element 
subregion,  the  order  of  reduction  and  the  number  of  elements  which  couple  the  external  region  with 
the  macro-element.  Since  this  number  is  much  bigger  in  3D  than  in  2D  domains,  the  performance 
of  the  reduction  algorithm  is  usually  significantly  affected  in  3D  problems.  In  order  to  overcome 
this  limitation  we  decrease  the  size  of  an  individual  macro-element  interface  by  projecting  the  fields 
of  its  boundary  onto  a  subspace  spanned  by  a  reduced  modal  base. 

The  modal  projection  and  the  reduction  process  used  in  our  approach  does  not  introduce  any 
frequency-dependent  elements  into  the  resultant  system  of  equations,  therefore  it  can  be  applied 
not  only  to  efficiently  solve  the  driven  problems  with  wideband  frequency  sweeps,  but  also  to 
compute  the  resonant  frequencies  of  the  cavities,  without  resorting  to  nonlinear  solvers. 
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Model  Order  Reduction  for  Systems  with  Geometrical  Parameters 
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With  Model  Order  Reduction  (MOR)  methods  an  original,  typically  large-scale  system  X  is  reduced 
to  a  smaller  system  X  which  approximates  the  transfer  function  of  X.  If  I]  is  parameter-dependent, 
it  is  wished  to  retain  the  dependence  in  £.  The  methods  for  the  latter  are  known  as  Paramet¬ 
ric  Model  Order  Reduction  (PMOR).  In  both  cases,  commonly  the  reduction  is  accomplished  by 
projection  (A.Antoulas,  ’Approximation  of  Large  Scale  Dynamical  Systems’,  SIAM,  2005). 

The  models  stem  from  the  Maxwell  grid  equations  which  are  obtained  from  the  continuous  Maxwell 
equations  with  the  help  of  the  Finite  Integration  Theory  (FIT)  (T.Weiland,  ’Time  domain  electro¬ 
magnetic  field  computation  with  finite  difference  methods’,  Int.J.Num.Mod,  9:295-319,1996).  For 
the  PMOR  described  in  the  following,  the  systems  are  derived  from  the  discrete  Helmholtz  equation 

£:  (M£(£)s2  +  C^itM-1(^)CFit)  x  =  Bi,  u  =  Cx,  (1) 

The  diagonal  matrices  M£,  contain  the  material  and  geometry  information  of  X  and  the  vector 
£  indicates  the  parameter  dependence.  The  topology-matrix  Cfit  represents  the  curl  operator. 
The  block- vectors  i  and  u  are  the  in-  and  output  of  X!  and  are  associated  to  the  matrices  B  and  C. 


The  most  important  challenge  in  calculating  appropriate  projection  matrices  V  are  geometrical 
variations,  as  the  models  usually  stem  from  automatically  created  meshes,  that  change  for  ge¬ 
ometry  variations.  Thus,  both  the  model  size  and  the  mesh  topology  change,  a  fact  that  most 
PMOR  methods  cannot  deal  with.  One  method  that  can  cope  with  geometrical  parameters  is  pre¬ 
sented  in  (K. Stavrakakis,  T. Wittig,  W. Ackermann,  T.Weiland,  ’Model  Order  Reduction  Methods 
for  Multivariate  Parameterized  Dynamical  Systems  obtained  by  the  Finite  Integration  Theory’, 
Proceedings  URSIGASS,  2011).  Projection  matrices  corresponding  to  local  systems  are  used 
to  form  Van  =  [Vi...Vn].  The  projection  matrix  V  is  chosen  as  the  leading  singular  vectors 
of  Vaii.  This  method  has  been  successfully  applied  to  parameterized  FIT  systems  (1),  as  shown 
in  the  above  reference.  In  order  to  work  on  a  constant  mesh  basis,  the  mesh  is  freezed  while  the 
parameter  is  varied.  Unlike  for  dielectric  material,  a  problem  arises  for  perfectly  electric  conducting 
(PEC)  material,  as  with  its  variation  the  zero  entries  in  the  material  matrices  change.  Typically, 
the  zeros  are  removed  in  order  to  avoid  redundancy.  Thus,  the  system  size  changes,  without  the 
mesh  being  changed.  This  can  be  handled  by  not  allowing  for  a  variation  greater  than  the  mesh-cell 
size.  Within  this  limit,  the  geometry  variation  is  carried  out  by  the  Partially  Filled  Cells  feature 
of  the  FIT.  An  interesting  case  which  needs  special  treatment  arises  if  the  variation  is  considered 
larger  than  the  mesh-cell  size. 
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Besides  frequency,  passive  microwave  structures  may  exhibit  a  number  of  design  variables,  such 
as  material  properties  or  geometric  parameters.  Many  practical  applications  require  the  solution 
response  of  such  structures  over  a  broad  frequency  band  and  a  wide  range  of  design  parameters. 

As  long  as  the  parameter-dependence  of  the  finite-element  (FE)  model  is  affine,  as  is  the  case 
for  operating  frequency  and  materials,  response  surfaces  are  constructed  very  efficiently  by  means 
of  projection-based  methods  of  model-order  reduction  (MOR):  robust  single-point  (SP)  and  multi¬ 
point  (MP)  algorithms  are  readily  available.  The  important  class  of  geometric  parameters,  however, 
does  not  lead  to  affine  parameter-dependence  and  is  hence  not  accessible  to  such  techniques.  To 
cope  with  this  situation,  several  methods  for  replacing  the  original  parameterization  by  some  affine 
interpolation  or  approximation  have  been  proposed,  e.g.  in  [1].  This  is  a  field  of  active  research, 
and  present  parametric  model-order  reduction  algorithms  [2],  [3]  exhibit  various  limitations,  e.g., 
with  respect  to  the  size  of  the  parameter  domain,  computational  complexity  or  convergence  rates. 

This  contribution  proposes  a  framework  for  generating  parametric  models  that  are  of  very  low 
dimension  and  thus  fast  to  evaluate.  It  is  capable  of  handling  geometric  design  parameters  as  well 
as  frequency  and  materials  and  features  convergence-rates  that  are  superior  to  competing  methods. 
In  addition,  important  system  properties,  such  as  passivity  or  reciprocity,  are  preserved.  The 
suggested  methodology  is  based  on  affine  parameter-reconstruction  combined  with  a  new  type  of 
parameter-dependent  projection  matrices. 

In  our  talk,  we  will  present  the  underlying  theory  and  give  numerical  examples  that  demonstrate 
the  characteristic  features  of  the  suggested  approach. 
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We  present  a  reduced  basis  element  method  (RBEM)  for  the  time-harmonic  Maxwell’s  equation 
after  introducing  the  reduced  basis  method  (RBM)  and  some  of  its  interesting  applications  including 
a  Pacman  scattering  problem.  RBM  is  indispensable  in  scenarios  where  a  large  number  of  numerical 
solutions  to  a  parametrized  partial  differential  equation  are  desired  in  a  fast /real-time  fashion. 
Thanks  to  an  offline-online  procedure  and  the  recognition  that  the  parameter-induced  solution 
manifolds  can  be  well  approximated  by  finite-dimensional  spaces,  RBM  can  improve  efficiency  by 
several  orders  of  magnitudes.  This  is  achieved  by  approximating  and  storing  the  solutions  of  the 
PDE  for  a  wisely-selected  set  of  parameters,  and  then  finding  the  reduced  basis  approximation 
for  a  new  parameter  as  a  Galerkin  projection  onto  the  low-dimensional  space  spanned  by  these 
pre-computed  solutions.  Its  accuracy  is  maintained  through  a  rigorous  a  posteriori  error  estimate. 

The  RBEM  is  RBM  with  a  particular  parameter,  that  is,  the  geometry  of  the  computational  domain, 
coupled  with  domain  decomposition  method.  The  basic  idea  is  to  first  decompose  the  computa¬ 
tional  domain  into  a  series  of  sub-domains  that  are  deformed  from  several  reference  domains;  then 
to  associate  with  each  reference  domain  pre-computed  solutions  to  the  same  governing  partial  dif¬ 
ferential  equation,  but  with  different  choices  of  deformations;  Finally  to  seek  the  approximation 
on  the  whole  domain  as  a  linear  combination  of  the  corresponding  pre-computed  solutions  on  each 
subdomain. 

We  will  introduce  a  new  RBEM  that  simulate  electromagnetic  wave  propagation  in  a  pipe  of  varying 
shape.  This  has  potential  application  in  fine-tuning  waveguide  design.  Unlike  the  traditional 
methods  with  domain  decomposition,  we  do  not  need  a  mortar  type  method  to  “glue”  the  various 
local  functions.  This  “gluing”  is  done  “automatically”  thanks  to  the  underlying  discontinuous 
Galerkin  method  we  are  using.  We  present  the  rationale  of  the  method  together  with  numerical 
results  showing  exponential  convergence. 
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We  present  reduced  basis  (RB)  approximation  and  associated  a  posteriori  error  estimation  for  rapid 
and  reliable  solution  of  parametrized  partial  differential  equations  (PDEs).  We  focus  in  particular 
on  saddle  point  problems  which  arise  in  a  wide  variety  of  applications  in  computational  science  and 
engineering.  Examples  of  sources  of  saddle  point  problems  include  mixed  finite  element  methods 
in  fluid  and  solid  mechanics,  and  optimization  with  PDEs.  This  particular  work  was  motivated 
by  and  applied  to  problems  in  incompressible  fluid  flow  through  parametrized  domains,  but  the 
methodology  can  be  applied  to  saddle  point  problems  in  general. 

Although  RB  methods  are  well-developed  for  several  classes  of  PDEs,  parametrized  saddle  point 
problems  pose  additional  difficulties  that  have  not  been  fully  addressed:  Parameter-dependent 
constraints  cause  complications  not  only  in  the  choice  of  stable  RB  approximation  spaces,  but  also 
in  the  construction  of  rigorous  and  computationally  efficient  a  posteriori  error  bounds.  Saddle 
point  systems  with  parameter-dependent  constraints  occur,  for  example,  when  applying  the  RB 
method  to  incompressible  fluid  flow  in  parametrized  domains.  Earlier  work  on  the  RB  method 
for  the  Stokes  and  incompressible  Navier-Stokes  equations  has  established  the  method  for  non- 
parametrized  domains  (i.e.,  for  the  case  of  parameter-independent  constraints).  There  have  also 
been  some  efforts  at  extending  the  method  to  the  Navier-Stokes  equations  in  parametrized  domains; 
however,  in  these  earlier  examples,  either  the  simple  geometric  variations  considered  are  applicable 
only  to  a  very  limited  set  of  problems,  or  rigorous  error  bounds  were  not  treated,  or  the  proposed 
error  bounds  require  the  evaluation  of  constants  which  are  very  difficult  to  compute. 

In  this  work,  we  focus  on  two  important  aspects.  First,  we  present  a  posteriori  error  estimates 
that,  unlike  earlier  approaches,  provide  separate  upper  bounds  for  the  errors  in  the  approximations 
for  the  primal  variable  and  the  Lagrange  multiplier  (e.g.,  the  velocity  and  pressure  variables,  re¬ 
spectively,  in  the  case  of  the  Stokes  equations).  The  proposed  method  is  a  direct  application  of 
Brezzi’s  theory  for  saddle  point  problems  to  the  RB  context,  and  exhibits  significant  advantages 
over  existing  a  posteriori  error  estimators  based  on  either  Babuskas  theory  for  noncoercive  problems 
or  a  penalty  approach.  Second,  in  view  of  approximation  stability  and  computational  cost,  we  shall 
analyze  and  compare  several  options  for  the  RB  approximation  spaces.  Through  numerical  tests, 
we  illustrate  the  significant  savings  (compared  to  earlier  strategies)  achievable  by  enriching  the  RB 
approximation  space  for  the  primal  variable  appropriately.  Finally,  both  a  posteriori  error  bounds 
and  enrichment  strategies  are  employed  in  an  adaptive  sampling  procedure  for  constructing  RB 
approximation  spaces  that  are  not  only  stable  but  also  efficient. 

We  present  numerical  results  for  a  problem  of  steady  Stokes  flow  through  a  two-dimensional  channel 
with  a  parametrized  rectangular  obstacle.  The  results  demonstrate  the  significant  effects  of  the 
enrichment  of  the  velocity  space  on  approximation  stability,  the  rapid  convergence  of  (stable) 
reduced  basis  approximations,  and  the  performance  of  the  proposed  error  bounds  with  respect  to 
sharpness  and  computational  cost. 
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Finite  Element  Methods  (FEM)  have  proven  to  be  robust  in  solving  the  time-harmonic  Maxwell’s 
equations  for  the  analysis  of  microwave  devices.  Nowadays,  microwave  engineering  relies  on  full- 
wave  analysis  not  only  to  accurately  predict  the  electrical  behaviour  prior  to  construction  of  mi¬ 
crowave  circuits,  but  also  to  actually  design  their  electrical  response.  In  this  last  framework, 
reduced-order  models  may  be  of  help,  since  many  modifications  will  be  carried  out  on  a  given  cir¬ 
cuit  either  by  experienced  engineers  or  global  optimization  techniques,  until  some  target  electrical 
response  is  achieved.  Indeed,  model  order  reduction  techniques  will  be  as  valuable  as  FEM  analysis 
if  the  reduced-order  model  is  a  reliable  surrogate  of  the  original  problem.  Otherwise,  we  may  not 
trust  on  these  techniques  for  design  automation,  although  they  are  computationally  more  efficient, 
and  may  still  rely  on  rather  time-consuming,  but  accurate,  FEM  analyses.  This  would  set  aside 
reduced-order  models  for  academic  purposes  but  never  for  industrial  use. 

In  this  work,  we  focus  on  a  rapid  and  reliable  frequency-parameter  sweep  in  microwave  circuits  based 
on  the  reduced-basis  method.  We  pay  special  attention  to  come  up  with  a  formulation  dropping 
singularities  in  frequency  out  of  the  Maxwell  operator.  The  proposed  model  order  reduction  is 
based  upon  the  following  observation:  should  the  field  solution  of  the  original  system  be  considered 
as  a  function  of  frequency,  the  field  solution  itself  does  not  arbitrarily  vary  as  frequency  changes, 
but  it  varies  smoothly  with  frequency,  i.e.,  the  field  solution  resides  on  a  very  low  dimensional 
subspace  induced  by  the  frequency  variation.  This  observation  is  a  direct  consequence  of  the  fact 
that  the  Maxwell  operator  as  well  as  the  forcing  term  are  smooth  functions  of  frequency.  The 
reduced-basis  method  focuses  on  approximating  the  evolution  of  the  electromagnetic  field  as  a 
function  of  frequency,  rather  than  approximating  the  electromagnetic  field  itself  as  the  FEM  does. 
We  propose  a  model  order  reduction  approach  where  field  solutions  at  different  frequencies  are  used 
as  new  basis  functions  to  capture  the  dynamics  of  the  electromagnetic  field  as  frequency  changes. 
Galerkin  projection  is  applied  to  reduce  the  original  FEM  system  into  a  small  system  detailing  the 
frequency  behavior  of  the  original  problem.  Thus,  the  frequency  response  of  the  structure  under 
analysis  can  be  obtained  with  ease.  Nevertheless,  we  still  need  to  solve  the  large  original  problem 
at  specific  frequencies  to  get  the  reduced  basis  functions.  In  this  work,  we  pay  special  attention 
to  adaptively  select  those  field  solution  frequency  samples  that  give  rise  to  a  robust  reduced  basis 
approximation  space.  In  addition,  the  accuracy  of  the  reduced-order  model  is  certified  by  means 
of  the  residual  error  of  the  reduced  field  solution  throughout  the  whole  frequency  band  of  interest. 
As  a  result,  a  fully  automatic  and  completely  reliable  model  order  reduction  process  is  achieved. 

Several  microwave  devices,  such  as  dieletric  resonator  filters,  dual-mode  filters  and  ultrawideband 
monopole  antennas,  will  illustrate  the  capabilities  of  this  approach. 
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The  Laboratoire  National  des  Champs  Magnetiques  Intenses  (LNCMI)  is  a  French  large  scale 
facility  enabling  researchers  to  perform  experiments  in  the  highest  possible  magnetic  field  (up  to 
35  T  static  field).  High  magnetic  fields  are  obtained  by  using  water  cooled  resistive  magnets  (cf. 


fig.  1)  connected  with  a  24  MW  power  supply  (C.  Trophime  et  ah,  Magnet  calculations  at  the 
GHMFL,  IEEE  Trans.  Applied  Superconductivity  12,  1,  1483-1487,  2002). 


Figure  1:  A  high  field  magnet:  (left)  exploded  view;  (center)  detailled  view  of  the  inner  part;  (right) 
temperature  of  the  inner  part  with  a  zoom  on  a  “turn” . 

The  design  and  optimization  of  these  magnets  requires  from  an  engineering  point  of  view  the  pre¬ 
diction  of  certain  “quantities  of  interest,”  or  performance  metrics,  which  we  shall  denote  outputs  — 
namely  magnetic  field  in  the  center,  maximum  stresses,  maximum  and  average  temperatures.  These 
outputs  are  expressed  as  functionals  of  field  variables  associated  with  a  set  of  coupled  parametrized 
partial  differential  equations  which  describe  the  physical  behavior  of  our  magnets.  The  parameters, 
which  we  shall  denote  inputs ,  represent  characterization  variables  —  such  as  physical  properties 
heat  transfer  coefficients,  water  temperature  and  flowrate,  and  geometric  variables  in  optimisation 
studies.  To  evaluate  these  implicit  input-output  relationships,  solutions  of  a  multi-physics  model 
involving  electro-thermal,  magnetostatics,  electro-thermal-mechanical  and  thermo-hydraulics  are 
requested.  It  should  be  noted  that  this  model  is  non-linear  as  the  material  properties  depend 
on  temperature.  In  practice  these  evaluations  represents  a  huge  of  computational  time  but  are 
mandatory  to  improve  the  magnet  design  as  we  cannot  rely  on  common  physical  sense. 

To  significantly  reduce  this  computational  time,  we  chose  to  use  model  order  reduction  strate¬ 
gies,  and  specifically  to  use  the  reduced  basis  method  which  is  well  adapted  to  the  evaluation 
of  input/output  relationships  (C.  Prud’homme,  D.  V.  Rovas,  K.  Veroy,  L.  Machiels,  Y.  Maday, 
A.  T.  Patera,  G.  Turinici,  Reliable  real-time  solution  of  parametrized  partial  differential  equa¬ 
tions:  Reduced-basis  output  bound  methods,  Journal  of  Fluids  Engineering,  124  ,  2002).  We 
will  present  the  reduced  basis  method  applied  to  the  non-linear  electro-thermal  coupled  prob¬ 
lem.  To  perform  this  work  we  rely  on  reduced  basis  framework  set  up  by  the  feel++  library 
(http:/ / www. feelpp.org/). 
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A  partially  orthonormal  high  order  basis  is  constructed  with  better  conditioning  properties  for 
tetrahedral  H(curl)-conforming  Nedelec  elements.  The  shape  functions  are  classified  into 
several  categories  with  respect  to  their  topological  entities  on  the  reference  3 -simplex.  The 
basis  functions  in  each  category  are  constructed  to  achieve  maximum  orthogonality.  The  new 
basis  the  composite  matrix  pM+S  has  better  conditioning  than  existing  basis. 
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Multi-physics  simulation  codes  are  designed  to  solve  systems  of  partial  differential  equations  (PDEs)  for 
modeling  phenomena  such  as  electromagnetic  wave  propagation,  shock-hydrodynamics,  elastic-plastic 
flow  and  magneto-hydrodynamics  (MHD).  In  this  talk  we  will  review  some  recent  work  at  LLNL  focused 
on  the  development,  software  implementation  and  application  of  high  order  finite  element  methods  for 
computational  physics.  The  common  theme  across  each  multi-physics  application  area  is  the  use  of  high 
order  finite  element  basis  functions  in  a  common  and  general  finite  element  discretization  approach. 

This  effort  started  in  the  context  of  the  linear  hyperbolic  Maxwell  equations,  where  it  is  known  that  high 
order  methods  are  capable  of  yielding  far  more  accurate  numerical  solutions  when  compared  to 
corresponding  grid  refined  2nd  order  methods,  and  often  at  a  significant  reduction  in  total  computational 
cost.  We  will  begin  with  a  review  of  the  development  of  high  order  H(Curl)  and  H(Div)  finite  element 
basis  functions  on  curvilinear  elements  combined  with  high  order,  energy  conserving  symplectic  time 
integration  algorithms  in  the  massively  parallel  CEM  code,  EMSolve.  We  will  review  their  use  in 
simulating  electrically  large  3D  wave  guiding  structures  and  highlight  some  benefits  of  this  approach. 
This  work  is  based  on  the  C++  software  library  FEMSTER  that  describes  an  object-oriented 
implementation  of  high  order,  curvilinear  finite  elements  at  the  element  level. 

Modifications  were  made  to  the  EMSolve  code  to  support  the  modeling  of  transient  magnetic  diffusion 
processes  under  the  so-called  MHD  approximation  of  Maxwell’s  equations.  We  will  review  the  use  of 
high  order  vector  finite  elements  for  solving  transient  magnetic  diffusion  problems  and  their  verification 
on  benchmark  /  analytic  problems.  This  work  was  ultimately  incorporated  into  existing,  massively 
parallel,  multi-material  arbitrary  Lagrangian-Eulerian  (ALE)  hydrodynamics  codes  for  the  purposes  of 
modeling  coupled  electro-thermal-mechanical  (or  MHD)  applications.  We  will  discuss  the  use  of 
divergence  preserving  finite  element  methods  in  the  context  of  magnetic  advection  and  diffusion  in  an 
ALE  setting  and  demonstrate  their  use  in  large-scale  3D  MHD  simulations. 

Finally,  we  consider  the  use  of  high  order  finite  element  methods  for  multi-material  Lagrangian 
hydrodynamics.  The  numerical  approximation  of  the  Euler  equations  of  gas  dynamics  in  a  moving 
Lagrangian  frame  is  at  the  heart  of  many  multi-physics  simulation  algorithms.  Recently,  we  have 
developed  a  general  framework  for  high-order  Lagrangian  discretizations  of  these  compressible  shock 
hydrodynamics  equations  using  high  order,  curvilinear  finite  elements.  This  work  is  based  on  our  current 
state  of  the  art,  open  source,  parallel  C++  library  for  finite  element  methods  research  and  development, 
called  MFEM  (mfem.googlecode.com).  We  conclude  the  talk  by  presenting  an  overview  of  this  approach 
and  its  implementation  in  the  BLAST  Lagrangian  hydro  research  code  (www.llnl.gov/casc/blast)  and  we 
demonstrate  the  advantages  that  high  order  curvilinear  finite  element  methods  can  provide  in  the  context 
of  (moving  mesh)  Lagrangian  hydrodynamics. 
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The  finite  element  method  (FEM)  is,  by  its  inherent  features,  especially  suitable  for  modeling  and 
analysis  of  structures  that  contain  inhomogeneous,  complex  electromagnetic  materials  and  geometrical 
irregularities.  The  FEM  is  well  established  as  a  method  of  choice  for  such  applications,  with  the  analysis 
of  open-region  scattering  structures  being  performed  truncating  the  FEM  domain  by  a  hybridization  with 
the  method  of  moments  (MoM)  or  by  some  sort  of  a  boundary  condition.  One  possible  general  strategy 
aimed  at  extending  the  practical  applicability  of  both  FEM  and  MoM  methods  over  their  inherent 
numerical  limits  and  considerably  enhancing  their  efficiency  in  real-world  simulations  is  the  diakoptic 
approach  (Olcan,  Ilic,  Notaros,  Kolundzija,  and  Djordjevic,  “Diakoptic  Higher-Order  FEM-MoM 
Approach,”  2010  IEEE  Antennas  and  Propagation  Society  International  Symposium  Digest),  according  to 
which,  the  solution  of  a  large  and  complex  electromagnetic  system  is  found  as  a  linear  combination  of 
solutions  of  diakoptic  subsystems,  using  explicit  linear  relations  between  coefficients  in  expansions  of 
equivalent  electric  and  magnetic  surface  currents  on  boundary  surfaces  of  subsystems. 

This  paper  presents  a  new  FEM-MoM-diakoptic  technique  for  scattering  analysis  in  the  frequency  domain 
that  is  specifically  designed  and  developed  for  higher  order  modeling  of  material  complexities  and  a  full 
exploitation  of  computational  efficiency  of  large  curved  finite  elements  with  p- refined  high-order  field 
approximations  in  applications  involving  arbitrary  material  anisotropy  and  inhomogeneity.  The  technique 
implements  Lagrange-type  generalized  curved  parametric  hexahedral  finite  elements  of  arbitrary 
geometrical-mapping  orders,  filled  with  anisotropic  inhomogeneous  materials  with  continuous  spatial 
variations  of  complex  relative  permittivity  and  permeability  tensors  described  by  Lagrange  interpolation 
polynomials  of  arbitrary  material-representation  orders.  Curl-conforming  hierarchical  polynomial  vector 
basis  functions  of  arbitrary  field-expansion  orders  are  used  for  the  approximation  of  the  electric  field 
vector  within  the  finite  elements,  while  fully  compatible  divergence-conforming  higher  order  vector  bases 
on  generalized  curved  parametric  quadrilaterals  are  implemented  for  both  MoM  and  diakoptic  surfaces. 
Being  solely  based  on  manipulations  with  expansion  coefficients,  the  diakoptic  approach  enables  great 
flexibility  in  selective  p-  and  /^-refinements  of  different  parts  of  a  complex  structure  combined  with 
entire-domain  modeling  of  anisotropic  inhomogeneous  electromagnetic  materials  -  for  the  best  overall 
performance  of  the  FEM-MoM-diakoptic  solver. 

The  FEM-MoM-diakoptics  is  validated,  evaluated,  and  discussed  in  a  variety  of  examples  of  scatterers, 
which  demonstrate  theoretical  features  and  computational  advancements  of  the  diakoptics  as  a  general 
modeling  methodology  and  in  conjunction  with  particular  higher  order  vector  basis  functions  and  curved 
elements  used  for  diakoptic  subsystems.  The  examples  show  benefits  of  using  truly  higher  order 
geometrical  and  current/field  modeling  of  both  electromagnetic  subsystems  and  diakoptic  surfaces,  when 
compared  to  low-order  basis  functions  and  elements.  They  also  show  the  efficiency  of  combining  the 
higher  order  bases  with  the  entire-domain  FEM  modeling  of  anisotropic  inhomogeneous  domains  with 
continuous  spatial  variations  of  material  parameters  (permittivity  and  permeability  tensors).  The  results 
demonstrate  that  the  more  complex  the  original  system  (problem),  in  terms  of  material  inhomogeneity  and 
anisotropy  and  geometrical  irregularity  of  diakoptic  subsystems,  the  more  advantageous  the  diakoptic 
analysis  over  the  standard  approaches  implementing  the  same  type  of  higher  order  numerical 
discretization. 
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The  FEniCS  Project  (www.fenicsproject.org)  is  a  collection  of  open  source  tools  focused  on  the 
automated  solution  of  differential  equations  using  the  finite  element  method  (FEM).  The  main  user 
interface  is  DOLFIN,  which  is  both  a  library  and  a  problem  solving  environment  that  provides  a 
simple  and  consistent  means  to  interact  with  the  other  FEniCS  components. 

One  of  the  main  strengths  of  FEniCS  is  that  it  provides  a  set  of  high-level  tools  with  which  to 
describe  a  finite  element  problem  and  ultimately  solve  the  resulting  equations.  This  paper  considers 
the  use  of  FEniCS  in  the  modelling  of  a  number  of  computational  electromagnetic  (CEM)  examples. 
This  is  aided  greatly  by  the  included  support  for  higher  order  basis  functions  from  Nedelec  and 
Lagrange  function  spaces  -  both  important  in  the  FEM  as  applied  to  CEM  -  that  allow  for  the 
rapid  implementation  and  testing  of  a  wide  range  of  formulation  aspects.  This  means  that  more 
time  can  be  spent  doing  new  research  as  opposed  to  concerning  oneself  over  implementation  issues. 

Although  FEniCS  provides  much  the  functionality  required  to  solve  basic  CEM  problems  such 
as  waveguide  cutoff  and  dispersion  analysis  (discussed  in  A.  Logg,  K.-A.  Mardal,  G.N.  Wells  et 
al.  ”  Automated  Solution  of  Differential  Equations  by  the  Finite  Element  Method,”  Volume  84 
of  Lecture  Notes  in  Computational  Science  and  Engineering,  Springer,  2011.),  as  well  as  resonant 
cavity  problems,  there  are  a  number  of  shortcomings  that  impede  its  use  for  more  advanced  prob¬ 
lems.  These  limitations  are  also  discussed  and  addressed  as  part  of  a  new  software  package  called 
SUCEM:FEM  (available  at  http://github.com/cemagg/sucem-fem)  which  is  built  on  FEniCS  and 
also  available  under  an  open  source  license.  This  new  solution  framework  allows  for  a  much  wider 
class  of  CEM  problems  -  including  radiation  and  scattering  problems  -  to  be  addressed. 
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It  is  well  known  that  in  the  vicinity  of  conducting  or  penetrable  edges  and  corners  certain  field 
components  and  the  surface  charge  and  current  densities  can  be  singular.  Although  this  behavior 
is  localized  at  the  edge  or  corner,  it  raises  the  cost  of  computation  since  one  needs  to  use  more 
expensive,  very  dense  meshes  in  the  neighborhood  of  the  singular  region.  The  literature  shows 
that  iterative  mesh-refinement  is  widely  used  in  the  FEM  context  to  treat  singularities,  which  is 
feasible  since  matrices  are  sparse  and  a  relatively  small  number  of  matrix  coefficients  must  be 
recomputed  for  the  refined  mesh.  However,  iterative  mesh  refinement  involves  complex  procedures 
and  codes,  uses  additional  unknowns  and  it  usually  results  in  an  increase  of  the  computational 
time  and/or  memory  requirements;  most  of  all,  any  iterative  mesh  refinement  technique  requires 
the  full  control  of  the  pre-processor  code(s)  used  to  discretize  the  computational  domain.  The 
best  alternative  for  the  accurate  treatment  of  fields  near  edges  is  the  introduction  of  singular  basis 
functions  of  the  additive  kind.  These  are  obtained  by  incrementing  the  regular  polynomial  vector 
bases  with  additional  degrees  of  freedom  to  model  the  singular  behavior  of  fields  and  currents. 
Interest  in  incorporating  edge  conditions  in  FEM  and  MoM  applications  dates  back  several  decades 
and  a  non- exhaustive  list  of  references  is  available  in  (Graglia  &  Lombardi,  IEEE  Trans .  Antennas 
Propagate  vol.  52,  pp.  1672-1685,  July  2004  and  Graglia  &  Lombardi,  IEEE  Trans.  Antennas 
Propagate  vol.  56,  pp.  3768-3788,  Dec.  2008).  This  presentation  will  review  the  development 
of  singular  basis  functions  for  applications  ranging  from  simple  scalar  problems  to  general  vector 
problems.  Some  of  the  issues  that  arise  include  the  following: 


•  The  number  of  degrees  of  freedom  associated  with  additive  basis  sets  formed  by  regular  and 
singular  basis  functions. 

•  Error  reduction  due  to  the  use  of  singular  basis  functions. 

•  Weight  functions  associated  with  singular  basis  function  subsets  and  implications  for  the 
commonly  used  Galerkin  testing  technique. 

•  Numerical  integration  techniques  used  to  compute  the  matrix  coefficients  arising  from  singular 
basis  functions. 

•  Metric  to  judge  the  relative  linear  independence  of  basis  sets  formed  by  regular  and  singular 
basis  functions. 

•  Deterioration  of  the  matrix  condition  number  due  to  the  singular  basis  function  subsets. 

•  Continuity  of  the  curl  or  divergence-conforming  singular  vector  functions  across  adjacent  cells. 

•  Treatment  of  corners  and  tips,  since  canonical  solutions  in  those  cases  are  as  not  well  known 
as  the  wedge. 
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It  is  of  critical  importance  to  efficiently  and  accurately  predict  global  resonances  of  a  3-D  electromagnetic 
structure  such  as  an  integrated  circuit  system  that  involves  arbitrarily  shaped  lossy  conductors  and 
inhomogeneous  materials.  A  quadratic  eigenvalue  solver  of  linear  complexity  and  electromagnetic 
accuracy  is  developed  in  this  work  to  fulfill  this  task.  Without  sacrificing  accuracy,  the  proposed 
eigenvalue  solver  has  shown  a  clear  advantage  over  state-of-the-art  eigenvalue  solvers  in  fast  CPU  time. 
It  successfully  solves  a  quadratic  eigenvalue  problem  of  over  2.5  million  unknowns  associated  with  a 
large-scale  3-D  on-chip  circuit  embedded  in  inhomogeneous  materials  in  40  minutes  on  a  single  3  GHz 
8222  SE  AMD  Opteron  processor. 

In  the  proposed  solver,  a  quadratic  eigenvalue  problem  is  first  formulated  for  analyzing  3-D  integrated 
circuit  structures  that  involve  arbitrarily  shaped  lossy  conductors  and  nonuniform  materials.  The  quadratic 
eigenvalue  problem  is  then  converted  to  a  generalized  eigenvalue  problem  to  facilitate  efficient 
computation.  The  computational  bottleneck  of  the  generalized  eigenvalue  problem  is  analyzed  and  found 
to  be  the  solution  of  a  large-scale  3-D  sparse  matrix.  A  set  of  orthogonal  prism  vector  bases  are  then 
developed  to  achieve  a  linear-complexity  solution  of  the  large-scale  sparse  matrix,  holding  the  complexity 
of  the  entire  eigenvalue  solution  to  linear.  Furthermore,  the  spurious  eigenvalues  associated  with  the 
orthogonal  vector  basis  based  discretization  of  the  generalized  eigenvalue  problem  are  identified.  The 
origin  of  the  spurious  eigenvalues  is  analyzed,  and  the  approach  to  removing  the  spurious  eigenvalues  is 
given.  In  addition,  the  accuracy  of  a  quadratic  eigenvalue  solution  is  investigated.  The  backward  error  is 
introduced  to  quantitatively  measure  the  accuracy  of  the  proposed  eigenvalue  solution.  An  optimal  scaling 
technique  is  used  to  transform  the  original  quadratic  eigenvalue  problem  to  a  scaled  quadratic  eigenvalue 
problem,  improving  the  accuracy  of  the  quadratic  eigenvalue  solution  by  a  few  orders  of  magnitude. 
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A  variety  of  physical  phenomena  such  as  acoustic  and  electromagnetic  time-harmonic  wave  propagation 
can  be  described  as  boundary-value  problems  governed  by  the  Helmholtz  equation.  It  is  well  known  that 
when  the  Helmholtz  operator  is  discretized  by  the  standard  Galerkin  finite  element  method  (FEM),  the 
accuracy  of  the  numerical  solution  deteriorates  rapidly  with  an  increasing  wavenumber,  and  the 
discretized  Helmholtz  operator  becomes  indefinite  (D.  Wang  et  al,  “Overview  of  the  discontinuous 
enrichment  method,  the  ultra- weak  variational  formulation,  and  the  partition  of  unity  method  for  acoustic 
scattering  in  the  medium  frequency  regime  and  performance  comparisons,”  Int.  J.  Numer.  Meth.  Engng ., 
2011,  in  press).  Normally,  a  higher  resolution  has  to  be  used  to  retain  a  certain  level  of  accuracy  at  the 
price  of  a  dramatic  increase  in  the  cost  of  the  standard  finite  element  analysis.  During  the  past  decades, 
several  FEMs  with  plane  wave  basis  functions  have  been  proposed  to  alleviate  the  negative  effect  of  the 
numerical  dispersion  to  achieve  an  efficient  finite  element  discretization  of  the  Helmholtz  equation.  These 
include  the  partition  of  unity  method  (PUM),  the  ultra-weak  variational  formulation  (UWVF),  and  the 
discontinuous  enrichment  method  (DEM)  (C.  Farhat,  I.  Harari,  and  L.  P.  Franca,  “The  discontinuous 
enrichment  method,”  Comput.  Methods  Appl.  Mech.  Engrg .,  vol.  190,  pp.  6455-4796,  2001).  Compared 
to  the  other  two  methods,  DEM  exhibited  higher  accuracy  and  better  computational  efficiency. 

Specifically,  the  enrichment  field  of  DEM  is  not  continuous  across  the  element  interfaces.  It  is  the 
Lagrange  multipliers  defined  on  element  boundaries  that  weakly  enforce  the  continuity.  Due  to  the 
discontinuity,  the  enrichment  can  be  eliminated  at  the  element  level  by  static  condensation,  which 
simplifies  the  method  and  improves  matrix  conditioning.  Furthermore,  if  combined  with  a  domain 
decomposition  method  (DDM),  the  above  condensed  algebraic  system,  which  is  related  to  element-level 
Lagrange  multipliers,  can  be  formulated  in  terms  of  the  subdomain  interface  Lagrange  multipliers  only. 
Similar  to  the  finite  element  tearing  and  interconnecting  (FETI)  method,  this  system  of  global  subdomain 
interface  Lagrange  multipliers  can  be  solved  iteratively  by  a  Krylov  subspace  method  (C.  Farhat,  R. 
Tezaur,  and  J.  Toivanen,  “A  domain  decomposition  method  for  discontinuous  Galerkin  discretizations  of 
Helmholtz  problems  with  plane  waves  and  Lagrange  multipliers,”  Int.  J.  Numer.  Meth.  Engng.,  vol.  78, 
pp.  1513-1531,  2009).  Therefore,  DEM  is  very  suitable  for  solving  large-scale  problems  on  parallel 
computing  facilities. 

One  important  research  topic  related  to  DEM  is  the  construction  of  plane  wave  basis  functions.  During 
the  past  few  years,  several  papers  have  been  published  on  the  choice  of  basis  functions  to  expand  the 
enrichment  field.  They  could  be  associated  with  either  quadrilateral/triangular  elements  for  two- 
dimensional  modeling,  or  hexahedral/tetrahedral  elements  for  three-dimensional  cases.  However,  all 
previous  work  is  limited  to  solving  the  scalar  two-  or  three-dimensional  Helmholtz  equation.  This  paper  is 
dedicated  to  the  implementation  of  DEM  for  solving  the  three-dimensional  vector  curl-curl  equation, 
which  is  the  governing  equation  for  electromagnetic  simulations.  The  potential  difficulties  include  the 
construction  of  vector  plane  wave  basis  functions,  which  are  the  nontrivial  solution  of  the  vector  curl-curl 
equation.  We  will  focus  on  basis  functions  compatible  with  triangular  and  tetrahedral  elements  since  they 
have  better  modeling  capability  for  arbitrarily  shaped  geometry.  We  will  demonstrate  the  efficiency  and 
scalability  of  the  discontinuous  enrichment  method  in  solving  three-dimensional  vector  wave  propagation 
and  scattering  problems. 
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Solving  wave  propagation  problems  with  FEM  results  in  a  large  number  of  unknowns  due  to  the 
huge  air  volume  to  be  modeled.  These  equation  systems  are  very  ill-conditioned  because  of  the  great 
material  differences,  element-size  changes  and  since  the  equation  systems  are  indefinite.  Common 
iterative  methods  (CG,  GMRES)  have  bad  convergence  quality  due  to  these  conditions. 

Applying  direct  solver  methods  to  overcome  the  problem  of  the  ill-conditioned  system  of  equations 
results  in  high  memory  requirements.  The  aim  of  this  paper  is  to  present  a  method  with  smaller  memory 
requirement  than  the  direct  methods  and  converging  faster  than  iterative  methods. 

A  possibility  of  decreasing  the  memory  requirement  of  the  common  LU  decomposition  is  using  a  Domain 
Decomposition  (DD)  based  Schur-complement  method.  Assembling  the  Schur  complement  matrix 
requires  a  long  time.  However,  solving  the  reduced  equation  system  doesn’t  need  the  full  Schur 
complement  matrix.  Applying  the  biconjugate  gradient  method  to  solve  the  Schur  complement  (DD- 
BiCG)  equation  system  will  result  in  an  efficient  iterative  solver  for  solving  wave  propagation  problems. 
The  efficiency  of  this  method  is  shown  on  a  dipole.  The  dipole  has  a  length  of  140  mm  a  width  of  1  mm 
and  thickness  of  20  pm.  There  is  a  160  pm  air  gap  in  the  middle.  An  air  volume  of  250  mm  is  modeled 
around  the  antenna.  The  air  volume  is  truncated  by  a  first  order  absorbing  boundary  condition  (ABC).  The 
voltage  is  prescribed  in  the  air  gap  (1  V,  1.5  GHz).  Modeling  an  eighth  of  the  problem,  using  A,v 
formulation  (A  is  the  magnetic  vector  potential,  v  is  the  modified  electric  scalar  potential),  and  second 
order  hexahedral  finite  elements  (20  nodes,  36  edges)  the  resulting  problem  has  1.986.152  edges  and 
669.398  nodes. 

The  efficiency  of  the  DD-BiCG  method  compared  with  the  incomplete  Cholesky  preconditioned 
Biconjugate  gradient  method  (IC-BiCG)  is  shown  in  Fig.  1.  The  iteration  will  be  finished  when  the  global 
relative  residual  is  smaller  than  10"7.  The  DD-BiCG  achieves  this  criterion  after  975  iterations  and  15.861 
seconds  using  41.8  GB  RAM.  The  IC-BiCG  was  not  able  to  converge  after  40.000  iterations  and  165.000 
seconds,  its  best  relative  residuum  is  1.198T0'5,  its  memory  requirement  is  5.25  GB  RAM.  With  the 
direct  method  PARDISO  implemented  in  Intel  FORTRAN  Math  Kernel  Library  required  81.5  GB  RAM 
and  hence  was  not  able  to  finish  the  calculation  on  a  64  GB  RAM  machine. 
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Fig.  1.  The  best  residuum  of  the  methods  during  the  iterations. 

Blue  line:  IC-BiCG;  Red  line:  DD-BiCG  residuum  on  the  interface;  Black  line:  DD-BiCG  residuum  in  the  whole  domain. 
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In  this  presentation  we  examine  in  detail  the  generation  of  quadrature  rules  for  singular  and  near-singular 
integrals  associated  with  integral  equations  and  higher  order  (curvilinear)  elements.  Though  the  approach 
is  geometry-representation  independent,  we  review  and  illustrate  it  using  the  Silvester-Lagrange 
interpolation  representation  of  curvilinear  triangles.  For  illustrative  purposes,  interpolator  divergence- 
conforming  bases  are  also  reviewed  and  assumed  (Graglia,  R.  D.,  D.  R.  Wilton,  and  A.  F.  Peterson, 
“Higher  Order  Interpolator  Vector  Bases  for  Computational  Electromagnetics,”  IEEE  Trans.  Antennas 
and  Propagate  Vol.  45,  Mar.  1997,  pp.  329-342). 

Given  an  observation  point  either  on  or  very  near  a  curvilinear  element,  one  should  first  construct  an 
appropriate  tangent  element.  The  point  of  tangency  (in  area  coordinates,  =  (^qj,^,^))  is  either 

the  observation  point  on  the  element,  or,  for  near  singularities,  the  closest  point  on  the  surface  or  its 
extension.  In  the  latter  case,  for  surfaces  convex  with  respect  to  the  observation  point,  the  closest  point  is 
found  by  solving  a  coupled  pair  of  non-linear  equations  via  Newton’s  method.  If  the  surface  is  concave 
with  respect  to  the  observation  point,  the  solution  may  not  be  unique. 

Position  vector  derivatives  with  respect  to  the  area  coordinates  determine  two  tangent  vectors,  A]  @  2  >  that 
also  form  two  edge  vectors  of  the  tangent  triangle;  the  third  is  merely ^3  —  —t\  —  ^2.  The  tangent 

triangle  vertex  locations,  element  Jacobian,  and  height  vectors  can  then  be  used  to  define  curvilinear  basis 
functions  at  the  tangent  point  in  terms  of  like  quantities  on  planar  triangles.  Also  points,  tangent  vectors, 
displacement  vectors,  and  Jacobians  on  the  tangent  triangle  differ  from  those  of  the  curvilinear  triangle 
only  to  second  order  with  respect  to  area  coordinates  about  .  Hence  singularity  cancellation  quadrature 

schemes  developed  for  handling  a  given  singularity  type  on  the  (planar)  tangent  triangle  are  valid  on  the 
curvilinear  triangle — only  a  few  more  sample  points  may  be  needed  to  account  for  the  higher  order 
variation  due  to  curvature.  Significantly,  for  divergence-conforming  bases,  the  Jacobians  appearing  in  the 
bases  and  differential  surface  elements  cancel,  so  that  all  curvature  effects  are  confined  to  the  edge 
vectors  ,  i  =  1, 2, 3,  and  the  displacement  term  |  r  —  r '  | . 
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Over  the  past  decade,  there  has  been  growing  interest  in  the  FDTD/FETD  hybrid  method  for  time- 
domain  electromagnetic  simulations.  This  method  models  part  of  the  computational  domain  where  extra 
modeling  accuracy  is  demanded  by  the  FETD  method,  while  applies  the  FDTD  method  elsewhere. 
Typically,  the  FDTD/FETD  hybrid  method  is  suitable  for  problems  with  locally  fine  features,  so  that  the 
computational  expensive  FETD  method  is  limited  to  certain  regions,  not  the  entire  domain. 

To  effectively  apply  the  FDTD/FETD  hybrid  method,  one  needs  to  generate  a  Cartesian/tetrahedral 
hybrid  mesh  by  embedding  the  tetrahedral  region  inside  background  Cartesian  grids.  That  is  to  say,  the 
boundary  of  the  tetrahedral  region,  which  is  the  FDTD/FETD  interface,  is  fixed.  Moreover,  the  tetrahedral 
region  is  desired  to  be  conformal  to  target  as  much  as  possible.  This  not  only  reduces  the  number  of 
degrees  of  freedom  (DOF)  in  the  FETD  analysis,  but  also  improves  the  computational  efficiency  of  direct 
matrix  solvers,  such  as  the  sparse  Cholesky  decomposition. 

Two  methods  are  common  in  automatic  mesh  generation  with  fixed  boundaries,  i.e.,  the  constrained 
Delaunay  and  the  advancing  front  technique.  It  was  observed  that  the  Delaunay  method  is  faster,  but 
meshes  generated  by  the  advancing  front  technique  are  often  of  higher  quality  and  much  easier  to 
improve.  Since  the  FETD  method  is  intended  to  be  applied  in  small  and  highly  confined  regions  in  the 
hybrid  approach,  mesh  quality  is  often  a  more  important  and  critical  issue  than  the  speed.  The  CPU  time 
of  mesh  generation  is  typically  much  less  than  that  of  hybrid  simulations,  while  mesh  quality  in  a  highly 
confined  region  is  difficult  to  achieve.  With  the  above  concerns,  the  advancing  front  technique  is  chosen 
for  hybrid  mesh  generation. 

In  this  presentation,  a  three-step  approach  was  presented  for  automatic  Cartesian/tetrahedral  hybrid 
mesh  generation.  With  a  prescribed  triangular  surface  model,  this  approach  "grows”  a  buffer  zone  from 
the  surface  of  a  target  at  first  and  positions  it  in  a  background  Cartesian  mesh.  The  tightness  of  the  buffer 
zone  can  be  adjusted  to  balance  the  needs  of  mesh  quality  and  computational  costs.  Subsequently,  a 
tetrahedral  mesh  is  generated  by  the  advancing  front  technique  with  respect  to  the  fixed  FDTD/FETD 
interface.  Finally,  the  quality  of  the  tetrahedral  mesh  is  further  improved  by  a  combination  of  edge 
splitting/merging,  face  and  edge  swapping,  smart  Laplacian  and  optimization-based  smoothing.  As  the 
result,  high-quality  meshes  with  good  computational  efficiency  in  the  subsequent  FDTD/FETD  hybrid 
analysis  can  be  generated. 
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The  finite  element  discretization  of  wave  equations  produces  a  numerical,  or  artificial,  dispersion  of 
propagating  waves.  As  the  Mixed  Finite  Element  Method  (MFEM)  is  based  on  the  discretization  of 
Maxwell’s  equations,  a  system  of  continuous  differential  operators,  it  retains  an  inherent  dispersive 
behavior.  The  objective  of  this  work  is  to  analyze  the  numerical  dispersion  of  the  mixed  FEM  for 
the  first-order  partial  differential  equations  based  on  Faraday’s  and  Ampere’s  laws  rather  than  the 
more  conventional  second-order  wave  equation. 

Two  different  approaches  are  applied  in  this  work  to  analyze  the  numerical  dispersion  of  MFEM. 
Both  approaches  assume  an  harmonic  plane  wave  propagating  through  an  infinite  two-dimensional 
vector  finite  element  mesh.  The  first  approach,  known  as  semi-discrete  analysis,  is  analytical,  and 
it  is  based  on  the  application  of  the  well-known  periodic  properties  of  an  harmonic  plane  wave 
to  a  portion  of  the  infinite  mesh;  thus,  it  obtains  an  expression  of  the  numerical  wave  number 
as  a  function  of  the  correct  value  and  direction  of  the  wavenumber.  This  approach  becomes  very 
cumbersome  for  higher  orders  basis  functions. 

A  second  approach,  known  as  modal  analysis,  requires  the  application  of  periodic  boundary  condi¬ 
tion  in  a  rectangular  domain  to  allow  propagation  of  plane  waves.  In  this  case,  a  plane  wave  with 
a  particular  wavenumber  (value  and  direction)  defines  an  approximate  eigenvector  for  the  system, 
hence  Rayleigh  quotient  is  exploited  to  obtain  an  approximate  numerical  wavenumber.  This  ap¬ 
proach  is  cross- validated  with  the  first  approach  for  low-order  basis  functions;  however,  high-order 
basis  functions  are  also  analyzed. 

Two  different  type  of  elements  are  analyzed:  triangular  and  quadrilateral.  For  triangular  elements, 
Mixed-order  Curl-Conforming  Edge-Based  Vector  basis  functions  (CtLn  and  LtQn)  are  used  (A. 
F.  Peterson,  S.  L.  Ray,  and  R.  Mittra,  Computational  methods  for  electromagnetics ,  IEEE  press 
New  York,  1998).  Quadrilateral  elements  adopt  spectral  basis  functions  based  on  Gauss-Lobatto- 
Legendre  polynomials,  with  an  arbitrary  order  of  interpolation  (J.  Lee,  T.  Xiao,  and  Q.  H.  Liu, 
“A  3-D  Spectral-Element  Method  Using  Mixed-Order  Curl-Conforming  Vector  Basis  Functions  for 
Electromagnetic  Fields”,  IEEE  Transactions  On  Microwave  Theory  And  Techniques,  Vol.  54,  No. 
1,  Jan  2006). 

To  understand  the  spurious  modes  produced  by  some  combinations  of  basis  function,  we  study 
the  effects  of  several  factors  on  the  numerical  dispersion:  direction  of  propagation,  order  of  basis 
functions,  type  of  element,  size  of  element,  and  variation  in  the  shape  of  element.  Finally,  the 
consequences  of  the  numerical  dispersion  on  some  parameters  in  frequency  and  time  domains  are 
also  analyzed,  such  as  phase  and  group  velocities,  phase  and  group  delays,  and  the  presence  of 
spurious  solutions  (non-monotonic  dispersion  curves). 
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To  engage  in  the  practical  design  of  electronic  components  and  systems,  one  must  contend  with  the 
presence  of  uncertainty  in  the  design  and  manufacturing  process  which  leads  to  material  and  structural 
variability  in  the  final  production  system.  To  maximize  production  yields  of  widely  multi-scale  structures 
under  conditions  of  uncertainty  in  a  high  volume  manufacturing  (HVM)  environment,  the  engineer  resorts 
to  electrical  design  automation  (EDA)  tools  in  the  search  for  a  robust  design  solution  space.  In  this 
context,  the  accuracy  and  efficiency  of  the  underlying  models  and  methods  of  EDA  are  fundamentally 
critical  to  the  characterization  efforts  of  the  component/system  and  in  providing  indispensable  guidance  to 
the  engineer  throughout  the  design  process.  To  enable  an  efficient  stochastic-based  design  optimization 
methodology  for  electrical  components  and  systems,  which  comprehends  random  material  and  geometric 
variability,  we  infuse  the  reduced-order  state-space  electromagnetic  model  together  with  the  stochastic 
collocation  method,  into  one  of  the  methods  used  extensively  for  electromagnetic  component/system 
analysis  and  design,  namely,  the  method  of  finite-difference  time-domain  (FDTD).  To  this  end,  we 
develop  the  stochastic  electromagnetic  macro-model  in  FDTD  by  formulating  an  abstraction  layer  that 
encapsulates  those  regions  where,  most  often,  uncertainty  is  present,  namely,  the  fine-features  of  the 
multi-scale  structure. 

In  essence,  the  electromagnetic  macro-model  is  a  reduced-order  state-space  representation  of  the 
discretized  Maxwell's  equations  for  the  FDTD  sub-domain  exhibiting  uncertainty.  From  the  full-order 
state-space  form,  it  is  possible  to  obtain  an  admittance  transfer  function  (ATF)  which  serves  to  map  the 
macro-model  and  the  FDTD  grid  at  the  union  of  their  common  boundary.  The  input  (excitation)  to  this 
ATF  is  the  magnetic  field  at  the  boundary  of  the  macro-model.  The  output  (response)  of  the  ATF  is  the 
electric  field  at  the  inner  perimeter  of  the  macro-model.  To  further  enhance  the  computational  efficiency 
of  operations  pertaining  to  the  system  ATF,  model  order  reduction  (MOR)  techniques  may  be  utilized  to 
reduce  the  degrees  of  freedom  of  the  full-order  state-space  model  while  maintaining  the  desired  accuracy. 
Thus,  we  apply  a  modified  version  of  the  enhanced  nodal  order  reduction  (ENOR)  technique  to  the  full- 
order  state-space  electromagnetic  system  ATF,  and  use  the  resulting  stochastic  projection  matrix  to  build 
the  passive  reduced-order  stochastic  electromagnetic  system.  We  name  this  reduced-order  stochastic 
abstraction  layer,  the  stochastic  electromagnetic  macro-model.  At  its  heart,  building  the  stochastic  macro¬ 
model  relies  on  highly  accurate  numerical  integration  over  random  space;  for  this,  we  utilize  the  sparse 
Smoljak  grid  and  the  Clenshaw-Curtis  quadrature  rule  in  one-dimensional  random  space,  and  extend  these 
principles  to  cubature  in  multi-dimensional  random  space. 

In  the  current  state  of  the  art  EDA,  and  in  absence  of  the  stochastic  electromagnetic  macro-model,  it  is 
necessary  to  perform  repeated  discretization  of  the  large  deterministic  domain  for  every  random  variation 
in  the  small  fine-featured  stochastic  sub-domains.  The  development  of  the  stochastic  electromagnetic 
macro-model  eliminates  the  need  for  repeated  discretization  of  the  overall  structure,  for  every  such 
random  variation.  Indeed,  only  a  single  FDTD  grid  needs  to  be  developed  for  the  deterministic  portion  of 
the  overall  structure  irrespective  of  the  realization  generated  by  a  specific  choice  of  the  random 
parameters  in  the  sub-domains  exhibiting  statistical  variability;  thus,  the  macro-model  results  in 
significant  computational  savings  by  eliminating  operations  pertaining  to  repeated  discretization  of  the 
deterministic  domain  for  each  variation  in  the  stochastic  sub-domains.  To  this  end,  we  provide  a 
methodology  for  the  efficient  generation  and  utilization  of  the  stochastic  macro-model  for  the  purpose  of 
time-domain  analysis  in  FDTD,  under  material/geometric  uncertainty  conditions. 
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The  use  of  time-domain  methods  to  simulate  electromagnetic  problems  is  well  rooted  in  the 
electromagnetics  community  due  to  several  reasons,  one  being  the  ability  to  obtain  wide  band  frequency 
results  in  a  single  simulation.  The  Finite-Difference  Time-Domain  (FDTD)  method  in  particular  is  one  of 
the  most  popular  methods  utilized  for  this  purpose.  However,  one  of  the  limiting  factors  inherent  to 
standard  FDTD  is  the  need  to  conform  a  grid  to  a  coordinate  system,  which  poses  significant  problems 
when  attempting  to  model  large  objects  with  fine  features  or  small  objects  in  large  domains.  One 
common  solution  proposed  to  overcome  this  difficulty  is  the  use  of  spatial  sub-gridding  techniques.  This 
solution  itself  poses  new  problems  such  as  instability  due  to  the  insertion  of  coarse/fine  grid  interface(s). 
Many  of  the  proposed  methods  for  spatial  subgridding,  while  handling  dissimilar  grid  interfaces 
effectively,  ignore  time-discretization  all  together.  This  imposes  further  restrictions  since  the  time-step 
size  allowable  in  order  to  maintain  numerical  stability  due  to  the  Courant  stability  limit  is  dictated  by  the 
most  finely  gridded  region.  This  dependence  on  a  global  time  step  is  sub-optimal  for  the  more  coarsely 
gridded  region(s)  causing  further  complications  such  as  numerical  dispersion. 

The  ultimate  goal  for  a  subgridding  scheme  would  be  for  each  dissimilarly  gridded  region  to  have  its  own 
optimal  time-step  size  regardless  of  the  granularity  of  the  remainder  of  the  domain.  The  method  proposed 
here  is  the  first  step  in  the  solution  to  address  this  issue  which  is  to  have  two  distinct  time-step  sizes  for  a 
uniformly  gridded  spatial  domain.  The  method  is  based  on  the  idea  that  a  regular  finite-element 
time-domain  (FETD)  hexagonal  mesh  can  be  proven  equivalent  to  standard  FDTD  method  through  the 
use  of  “lumping”  of  the  mass  matrices  or  simply,  mass-lumping.  The  purpose  for  deriving  an  FDTD 
scheme  via  finite-elements  is  an  effort  to  maintain  proper  field  continuities  and  physical  field  properties. 
The  procedure  is  then  extended  to  include  the  time  axis  and  is  considered  as  a  separate  and  discretizable 
axis  similar  to  the  spatial  axes  and  is  “time-lumped”  accordingly.  The  lumped  FETD  elements  are  then 
assembled  via  assembly  or  joining  matrices  and  FDTD  stencils  with  special  updates  along  the  time 
boundary  are  derived  and  introduced  into  a  standard  FDTD  scheme. 


The  method  described  above  was  applied  to  a  two-dimensional  FDTD  grid  spatially  discretized  uniformly 
throughout  the  domain  and  divided  into  two  regions.  One  portion  of  the  grid  was  time-stepped  at  one 
temporal  discretization  size  (coarse-time  grid),  while  the  other  portion  of  the  grid  (fine-time  grid)  is 
stepped  at  a  rate  double  of  the 
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coarse-time  step  size.  The  domain 
boundaries  were  terminated  with 
perfect  electric  conductor  (PEC) 
material.  As  a  preliminary 
indicator  for  stability  the  domain 
was  excited  with  a  Gaussian  pulse 
which  was  allowed  to  propagate 
for  over  300,000  time-steps.  At 
every  time  step  the  discrete  energy 
in  the  electric  field  throughout  the  FDTD  grid  was  calculated  and  recorded.  As  the  figure  above  shows, 
for  the  proper  choice  of  Courant  limit,  the  mean  discrete  energy  did  not  increase.  This  suggests  that  the 
scheme  is  conditionally  stable  and  conservative. 
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P-adaption  -  increasing  the  polynomial  orders  of  elements  automatically  based  on  estimated  errors  -  is  in 
many  respects  preferable  to  the  alternative  of  refining  the  finite  element  mesh  (h- adaption).  It  avoids  the 
complexities  of  mesh  refinement  and,  when  the  fields  are  smooth,  provides  excellent  convergence.  By 
using  singular  elements,  it  can  even  perform  well  when  the  field  is  singular,  e.g.,  at  sharp  metal  edges. 
However,  it  is  ultimately  limited  by  the  maximum  available  polynomial  order.  Combined  hp- adaption 
offers  the  possibility  of  even  faster  convergence,  with  no  built-in  limit.  This  paper  explores  the 
application  of  hp- adaption  to  the  computation  of  the  scattering  parameters  of  microwave  devices  using 
hierarchical,  tetrahedral  finite  elements. 

In  true  /^-adaption,  after  every  solution  the  mesh  may  be  refined  in  some  places  and  the  polynomial  order 
increased  in  other  places.  This  requires  that  there  is  some  way  of  tracking  the  polynomial  order 
everywhere  in  the  mesh  through  successive  mesh  refinements.  For  example,  if  at  some  stage  the  order  of 
an  element  is  2  and  then  the  mesh  in  that  region  is  refined,  the  order  2  has  to  be  carried  to  some  of  the 
new  elements  created.  This  presents  a  problem  for  many  mesh  refinement  schemes,  in  which  an  element 
created  by  refinement  does  not  necessarily  fit  inside  the  boundaries  of  an  earlier  coarse  element  -  what 
order,  then,  should  be  assigned  to  the  new  element? 

One  way  around  this  problem  is  to  use  special  elements  with  hanging  nodes.  These  permit  a  form  of  h- 
refinement  is  which  each  new  element  “belongs”  to  a  coarser  parent,  so  there  is  no  ambiguity  about  the 
order  of  the  refined  element.  A  good  example  is  (L.  E.  Garcia-Castillo,  D.  Pardo,  L.F.  Demkowicz, 
’'Energy-norm-based  and  goal-oriented  automatic  hp  adaptivity  for  electromagnetics:  application  to 
waveguide  discontinuities,"  IEEE  Trans.  MET ,  vol.56,  no.  12,  pp. 3039-3049,  Dec.  2008).  An  alternative  is 
proposed  here.  Any  set  of  hierarchical,  tetrahedral  elements  can  be  used,  but  a  mesh  refinement  scheme  is 
employed  that  preserves,  to  a  large  extent,  the  nested  character  that  is  needed  for  /^-adaption.  The 
scheme  is  due  to  J.  Bey  ("Tetrahedral  grid  refinement,"  Computing ,  vol.  55,  pp.  355-378,  1995).  The 
essential  refinement  operation  is  the  subdivision  of  a  tetrahedron  into  8  “sons”.  In  addition,  there  are 
several  “irregular”  refinements  that  permit  a  transition  from  a  coarse  region  to  a  more  finely  subdivided 
region.  The  sons  produced  by  irregular  refinement  are  temporary:  when  one  of  them  needs  to  be  refined, 
all  of  them  are  removed  and  the  parent  tetrahedron  is  refined  regularly.  In  this  way,  there  is  no 
propagation  of  irregular  shapes. 

An  hp- adaptive  scheme  based  on  Bey  mesh  refinement  is  applied  to  the  computation  of  scattering 
parameters  in  the  frequency  domain.  An  error  indicator  previously  used  for  /^-adaption,  which  targets  the 
error  in  the  scattering  parameters,  is  used  to  guide  the  adaption  (D.  Nair,  J.  P.  Webb,  "P-adaptive 
computation  of  the  scattering  parameters  of  3-D  microwave  devices,"  IEEE  Trans.  Magnetics ,  vol.  40,  pp. 
1428-1431,  Mar  2004). 
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Efficient  RCS  prediction  of  jet  engine  air  intakes  using  /ip-FEM 

and  reduced  order  modelling 
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Jet  engine  inlets  are  often  electrically  huge  semi-open  channels.  Measurements  have  revealed  that 
these  channels  are  one  of  the  major  contributors  to  the  overall  RCS  signature  of  an  airplane. 
Reliable  and  accurate  Radar  Cross-Section  (RCS)  prediction  of  air  intakes  remain  one  of  the  most 
challenging  problems  in  computational  electromagnetics,  (CEM).  For  low  RCS  they  are  often  curved 
and  partially  coated  with  radar  absorbing  materials  (RAM)  to  hide  jet  engine  fan  disks  at  the  far 
end,  see  Figure  1. 

Highly  efficient  parallelised  computational  techniques  have  to  be  used.  Multi-level  domain- 
decomposition  (DD)  techniques  provide  an  efficient  way  to  divide  and  conquer  the  resulting  huge 
computational  problem.  Different  models  may  be  used  and  mixed  along  the  channel  to  obtain  the 
requested  efficiency,  accuracy  and  reliability  of  the  RCS  prediction.  We  use  the  so-called  hp-ve rsion 
of  the  Finite  Element  Method  (FEM)  for  modelling  sections  with  material  inhomogeneities  and  for 
modelling  the  scattering  on  the  terminating  section  consisting  of  disks  with  irregular  jet  engine  fan 
blades.  We  discuss  cavity  RCS  predictions  based  on  a  reduced  order  uncoupled  interior  problem 
analysis  and  the  Kirchhoff  aperture  integration  procedure.  The  efficiency  gain  and  modelling  error 
introduced  using  this  simplified  approach  is  discussed.  The  DD-technique  in  terms  of  Neumann-to- 
Dirichlet  (NtD)  operators  is  outlined.  The  sub-domain  wise  contributions  are  assembled  projecting 
the  contributions  of  coupling  neighbours  onto  an  auxiliary  complete  basis  at  the  common  interface. 
A  choice  we  investigate  for  a  homogeneous  channel  cross-section  is  the  (TE-,TM-)  waveguide  modal 
basis.  The  use  of  truncated  projections  is  investigated. 

The  efficiency  of  the  method  is  illustrated.  The  reduced  order  hp- FE  sub-domain  based  scattering 
matrix  methodology  is  proven  suitable  for  calculating  the  Radar  Cross-Section  (RCS)  for  electrically 
huge  jet  engine  air  intakes.  The  efficiency  gain  in  degrees  of  freedom  obtained  by  using  hp-ve  rsion 
FEM  instead  of  classical  low  order  h- version  FEM  is  shown  to  be  roughly  one  order  of  magnitude. 
The  model  reduction  achieved  by  changing  from  inter-facial  FE-d.o.f:s  to  guided  wave  participation 
factors  implies  another  gain  in  degrees  of  freedom  which  becomes  very  substantial  for  air  intakes 
with  electrically  large  homogeneous  sections.  It  is  shown  that  the  modal  reduction  can  be  made 
without  significant  loss  of  accuracy  in  the  cavity- RCS. 


Figure  1:  Low  RCS  S-shaped  jet  engine  air  intake. 
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There  has  lately  been  a  noticeable  interest  within  the  computational  electromagnetics  (CEM)  community 
for  higher  order  techniques,  which  rely  on  basis  functions  of  high  orders  set  on  electrically  large 
(whenever  possible)  elements  (B.  M.  Notaros,  “Higher  Order  Frequency-Domain  Computational 
Electromagnetics,”  IEEE  Trans.  Antennas  Propagat.,  Vol.  56,  August  2008).  Additionally,  higher  order 
techniques  allow  for  efficient  utilization  of  irregular,  elongated,  and  curved  elements,  defined  using 
Lagrange  interpolating  polynomials,  Bezier  curves,  or  NURBS  curves,  which  can  be  tremendously 
beneficial  in  creating  effective  meshes  and  in  yielding  the  final  systems  of  equations  with  significantly 
fewer  number  of  unknowns.  Most  importantly,  higher  order  bases  and  elements  enable  efficient  p-  and 
/^-refinements,  as  well  as  adaptive  CEM  schemes. 

However,  the  great  modeling  flexibility  of  higher  order  elements,  basis  and  testing  functions,  and 
integration  procedures,  which  is  the  principal  advantage  of  the  higher  order  CEM,  is  also  its  greatest 
shortcoming.  Namely,  it  poses  numerous  dilemmas,  uncertainties,  options,  and  decisions  to  be  made  on 
how  to  actually  use  these  elements  and  functions.  In  other  words,  with  the  additional  degrees  of  freedom 
in  modeling,  a  user  has  to  handle  many  more  parameters  in  building  a  CEM  model,  which  requires  a 
great  deal  of  modeling  experience  and  expertise,  and  possibly  considerably  increases  the  overall 
simulation  (modeling  plus  computation)  time.  Examples  of  questions  that  need  to  be  addressed  are  how 
large  the  elements  can  be,  what  polynomial  orders  of  bases  should  be  used  in  particular  cases,  and  how 
accurate  the  integration  has  to  be  depending  on  the  adopted  polynomial  orders  in  different  directions. 

Within  efforts  to  answer  these  questions,  we  have  recently  established  and  validated  general  guidelines 
and  instructions,  and  as  precise  as  possible  quantitative  rules,  for  adoptions  of  optimal  higher  order 
parameters  for  electromagnetic  modeling  using  the  method  of  moments  (MoM)  in  single  precision 
computations.  For  instance,  we  have  concluded  that  (in  cases  with  well-behaved  fields)  the  Lagrange-type 
(curved  or  fiat)  quadrilaterals  as  large  as  up  to  two  wavelengths  on  a  side,  with  polynomial  orders  of  the 
equivalent  surface  current  approximation  of  N  =  6  in  each  direction  and  the  number  of  points  in  Gauss- 
Legendre  quadrature  integrations  of  about  N  +  2,  are  generally  optimal. 

In  our  continued  study  of  higher  order  parameters  in  CEM,  this  paper  addresses  the  same  problem  and 
defines  similar  rules  for  the  finite  element  method  (FEM)  using  p-  and  @-refined  higher  order 
hierarchical  generalized  hexahedral  elements  and  models  in  double  precision  computations.  To  define 
these  rules  and  draw  some  general  conclusions,  we  carefully  and  systematically  study  a  diverse  and 
comprehensive  set  of  higher  order  FEM  simulations  of  three-dimensional  (3-D)  cavities,  2-D  waveguides, 
and  3-D  scatterers  (with  the  FEM  domain  closed  using  the  first  order  absorbing  boundary  condition  or  a 
fixed  high-order  boundary  integral),  applying  exhaustive  sweeps  in  frequency,  polynomial  orders  of  the 
field  approximation,  and  numbers  of  points  in  the  Gauss-Legendre  integration  (integration  accuracy). 
Preliminary  results  show  that,  while  it  can  generally  be  adopted  that  curved  Lagrange-type  hexahedral 
finite  elements  can  be  as  large  as  two  wavelengths  on  a  side  (analogously  to  the  MoM  analysis), 
somewhat  larger  numbers  of  integration  points  can  be  used  in  each  direction  than  those  reported  in  the 
MoM  study.  However,  higher  integration  accuracy  results  in  possibly  substantially  longer  matrix  filling 
times;  hence,  in  this  work,  we  also  discuss  a  compromise  between  the  two  requirements. 
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When  high  order,  hierarchical,  tetrahedral  finite  elements  are  used  to  solve  large,  time-harmonic, 
electromagnetic  wave  problems,  an  efficient  way  to  solve  the  resulting  matrix  equation  is  to  use  a  Krylov 
iterative  method  with  a  multilevel  preconditioner,  e.g.,  P.  Ingelstrom,  V.  Hill,  R.  Dyczij-Edlinger, 
"Comparison  of  hierarchical  basis  functions  for  efficient  multilevel  solvers,"  IET  Science  Measurement  & 
Technology ,  vol.  1,  pp.  48-52,  Jan  2007.  The  preconditioner  uses  a  few  steps  of  Gauss-Seidel  relaxation  at 
the  second  order  and  higher  levels,  but  requires,  at  each  iteration,  the  direct  solution  of  a  matrix  problem 
corresponding  to  the  lowest  order  level,  i.e.,  Whitney  elements  with  one  degree  of  freedom  per  edge  of 
the  mesh.  For  large  problems  with  fine  geometric  detail,  the  number  of  elements  may  be  large  enough  that 
direct  solution  on  this  level  requires  too  much  memory  or  takes  too  long.  An  iterative  alternative  is 
desirable.  One  possibility  that  has  been  widely  explored  is  domain  decomposition.  We  consider  an 
alternative  method  that  does  not  require  domain  decomposition. 

The  simplest  approach  is  to  replace  the  direct  solution  by  a  single  step  of  Gauss-Seidel  relaxation. 
Actually,  to  take  into  account  the  kernel  of  the  curl  operator  it  is  necessary  to  use  what  is  known  as  a 
hybrid  smoother,  in  which  a  second  relaxation  is  applied  in  a  space  corresponding  to  nodal  scalar  basis 
functions.  This  approach,  PI,  does  work,  but  the  number  of  iterations  is  considerably  greater  than  with 
direct  solution. 

An  alternative  (P2)  is  to  apply  an  auxiliary  space  preconditioner  (R.  Hiptmair,  J.  C.  Xu,  "Auxiliary  space 
preconditioning  for  edge  elements,"  IEEE  Transactions  on  Magnetics ,  vol.  44,  pp.  938-941,  Jun  2008). 
To  date,  this  has  been  applied  only  to  quasti-static  problems  and  has  not  been  used  with  higher  order 
elements.  It  requires  the  solution  of  two  auxiliary  nodal  problems:  the  first  one  corresponding  to  nodal 
scalar  functions,  as  in  PI,  and  the  second  corresponding  to  nodal  vector  functions.  The  nodal  problems 
are  solved  by  an  algebraic  multigrid  method.  We  employ  a  W-cycle,  solving  each  nodal  problem  three 
times  per  Krylov  iteration. 

The  table  below  shows  the  results  for  a  sample  problem,  an  E-plane,  corrugated,  waveguide  filter,  with  an 
unstructured  mesh  of  226,160  tetrahedra.  The  waveguide  cross  section  is  19.0mm  x  9.5mm  and  the  filter 
length  is  91.5mm,  which  is  about  3  wavelengths  at  the  analysis  frequency  (10  GHz).  To  solve  each  nodal 
problem,  three  levels  of  algebraic  multigrid  are  used,  followed  by  a  direct  solve  at  the  coarsest  level 
(matrix  dimension  5,665).  The  new  preconditioner  (P2)  substantially  reduces  the  number  of  iterations. 
Notice  that  at  order  3  the  reduction  in  the  number  of  iterations  has  more  impact  because  of  the  extra  cost 
associated  with  the  second  and  third  level  parts  of  the  overall  preconditioner  and  this  is  reflected  in  the 
greater  speed-up. 


Element 

order 

Matrix 

dimension 

Preconditioner 

Krylov 

iterations 

CPU  time 
(seconds) 

Speed-up 

i 

352,971 

pi 

370 

106 

i 

352,971 

P2 

54 

47 

2.3 

3 

4,415,571 

PI 

349 

2,790 

3 

4,415,571 

P2 

60 

531 

5.2 
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This  presentation  will  continue  an  exploration  of  limitations  and  possibilities  for  discrete  differential 
forms,  and  in  particular,  the  Whiteny  forms  which  generalize  edge  and  face  elements  which  made  their 
way  into  the  finite  element  analysis  of  electromagnetic  fields  in  the  1980s.  Starting  in  the  1950s,  Andre 
Weil,  H.  Whitney,  H.  Federer,  D.  Sullivan,  J.  Dodzuik,  J.  Dupont,  W.  Muller,  R.  Hiptmair,  S.  Wilson,  E. 
Getzler  and  P.  Mnev,  contributed  to  the  development  and  application  of  Whitney  forms  to  various 
problems  in  topology  and  computational  electromagnetics  and  more  recently  to  computational  fluid 
dynamics  (the  most  recent  work  can  be  found  by  searching  on  http://arXiv.org  ). 

This  presentation  will  build  on  the  author’s  presentation  of  the  de  Rham  and  Whitney  maps 
(COMPUMAG  2007;  IEEE  Trans  Mag,  2008)  in  order  to  characterize  the  kernel  of  the  Whitney  map  and 
explore  the  implications  for  computational  electromagnetics.  Specifically,  the  contributions  of  the  above 
authors  are  considered  in  order  to: 

1)  Standardize  notation  between  application  areas  in  order  to  facilitate  an  exchange  of  results. 

2)  Develop  a  framework  for  comparing  error  analyses  in  various  fields.  For  example,  the  Sobolev 
space  framework  in  finite  element  analysis,  the  norms  of  geometric  measure  theory  (GMT)  and 
the  IR-UV  decompositions  of  quantum  field  theory  will  be  contrasted. 

3)  Develop  frameworks  for  both  h-p  refinement  in  the  finite  element  analysis  of  electromagnetic 
fields  and  measures  of  simplicity  of  “cuts”  for  magnetic  scalar  potentials. 

In  all  these  aspects,  the  Dupont  chain  homotopy  operator  plays  a  central  role.  How  this  operator 
characterizes  the  kernel  of  the  de  Rham  map  and  how  it  is  underappreciated  in  the  computational 
electromagnetics  community  will  be  stressed. 
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The  aim  of  this  work  is  to  compare  the  computational  performance  of  nodal-based  finite  element 
formulations  versus  edge-based  finite  element  formulations.  Edge-based  formulations  are  widely 
adopted  in  the  EM  community  and  known  to  be  stable  (M.  Salazar-Palma  et  ah,  ’’Iterative  and 
Self-Adaptive  Finite-Elements  in  Electromagnetic  Modeling”.  Artech  House  Publishers,  1998.  J. 
M.  Jin,  ’’The  Finite  Element  Method  in  Electromagnetics”,  John  Wiley  &  Sons,  2002).  However, 
straightforward  nodal-based  finite  element  formulations  suffer  from  stability  problems.  Recently, 
two  nodal  based  formulations  overcoming  the  mentioned  stability  problems  have  been  proposed 
by  the  authors.  Details  of  the  two  formulations  are  given  below.  The  comparison  is  established 
between  each  one  of  the  nodal-based  formulations  and  their  edge-based  counterparts. 

The  first  nodal  approximation  is  the  simplification  of  the  weighted  regularized  Maxwell  equation 
method  explained  in  (R.  Otin.  ’’Regularized  Maxwell  equations  and  nodal  finite  elements  for  elec¬ 
tromagnetic  field  computations.”  Electromagnetics,  vol.  30,  pp.  190-204,  2010).  This  formulation 
presents  several  advantages:  It  provides  spurious-free  solutions  and  well-conditioned  matrices  with¬ 
out  the  need  of  Lagrange  multipliers,  its  integral  representation  is  well-suited  for  hybridization  with 
integral  numerical  techniques  and  the  nodal  solution  of  the  electromagnetic  domain  is  usually  easier 
to  couple  in  multi-physics  problems.  On  the  other  hand,  this  regularized  formulation  presents  a 
serious  drawback:  if  the  electromagnetic  field  has  a  singularity  in  the  problem  domain,  a  globally 
wrong  solution  is  obtained.  Fortunately,  this  disadvantage  can  be  corrected,  but,  in  doing  this, 
the  well-conditioning  of  the  matrix  can  be  put  in  jeopardy.  Then,  the  question  we  try  to  solve 
in  this  work  is  if  the  computational  advantages  of  the  proposed  nodal  formulation  compensate  its 
drawbacks  when  compared  with  the  well-establish  edge  double-curl  finite  element  formulation. 

The  second  nodal  formulation  we  are  going  to  test  is  based  on  the  algorithm  explained  in  (S.  Ba¬ 
dia  and  R.  Codina,  ”A  nodal-based  finite  element  approximation  of  the  Maxwell  problem  suitable 
for  singular  solutions,”  submitted )  for  the  Maxwell  operator.  This  formulation  will  be  adapted 
to  the  frequency  domain  and  compared  with  mixed  edge  element  formulations.  The  main  differ¬ 
ences  between  this  new  nodal  formulation  and  the  one  based  on  weighted  regularization  is  the 
fact  that  a  Lagrange  multiplier  is  needed  in  order  to  enforce  the  divergence-free  constraint  at  the 
discrete  level  but  there  is  no  need  to  define  the  weighting  functions.  This  last  point  is  important, 
since  it  does  not  require  to  determine  a  priori  were  the  singularities  are,  making  the  method  eas¬ 
ier  to  apply  and  implement.  On  the  other  hand,  it  does  not  exhibit  the  problem  of  integrating 
numerically  the  weighted  terms  and  no  additional  (regularized)  boundary  conditions  are  needed. 
Finally,  this  method  has  been  easily  extended  to  coefficient  jumps  in  (S.  Badia  and  R.  Codina, 
”A  combined  nodal  continuous-discontinuous  finite  element  formulation  for  the  Maxwell  problem,” 
Applied  Mathematics  and  Computation,  vol.  218,  pp.  4276-4294,  2011)  using  dG-type  terms  on 
the  material  interfaces. 
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The  Finite  Element  Method  (FEM)  enables  the  use  of  ’’adapted”  meshes,  not  only  to  the  geometry 
of  the  problem  domain,  but  to  the  solution  of  the  problem  itself.  One  of  the  most  powerful  methods 
to  generate  these  ’’adapted”  meshes  is  through  /ip-adaptivity  (i.e.,  by  simultaneously  varying  the 
size  h  and  the  polynomial  order  p  of  the  finite  elements  of  the  mesh),  hp- adaptivity  provides 
exponential  rates  of  solution  convergence,  even  in  the  presence  of  singularities,  in  contrast  to  h 
and  p  schemes,  in  which  algebraic  rates  of  solution  convergence  are,  in  general,  obtained.  Thus, 
very  accurate  solutions  can  be  obtained  with  an  hp-adaptive  strategy,  even  in  the  presence  of 
singularities.  Equivalently,  approximate  solutions  within  engineering  accuracy  can  be  obtained 
using  a  minimum  number  of  unknowns. 

In  the  last  years,  the  authors  have  developed  two  dimensional  (2D)  automatic  hp-adaptive  FE 
methods  for  electrodynamic  (i.e.,  full  wave)  problems.  The  problems  solved  include  closed  domain 
and  open  region  problems;  specifically,  waveguiding  problems  (L.  E.  Garcia-Castillo  et  ah,  IEEE 
Trans.  MTT,  12,  3039-3049,  2008),  radiation  and  RCS  computation  problems  (I.  Gomez-Revuelto 
at  ah,  IEEE  Trans.  MAG,  43(4),  1337-1340,  2007),  and  geophysics  applications  (D.  Pardo  et  ah, 
Comp.  Meth.  App.  Mech.  Eng.,  197,  3836-3849,  2008). 

In  this  paper,  a  three-dimensional  (3D)  fully  automatic  energy-norm  based  Ap-adaptive  Finite  Ele¬ 
ment  (FE)  strategy  is  described  and  applied  to  solve  relevant  electromagnetic  waveguide  structures. 
It  is  based  on  the  work  on  hp  adaptivity  leaded  by  Prof.  Demkowicz  of  the  University  of  Texas  at 
Austin.  As  in  the  2D  case,  the  three  dimensional  hp-strategy  supports  anisotropic  refinements  on 
irregular  meshes  with  hanging  nodes,  and  isoparametric  as  well  as  exact-geometry  elements. 

Extensive  numerical  results  demonstrate  the  suitability  of  the  Ap-adaptive  method  for  solving  dif¬ 
ferent  rectangular  waveguide  discontinuities.  These  results  illustrate  the  flexibility,  reliability,  and 
high-accuracy  of  the  method.  Solutions  of  waveguide  problems  delivered  by  the  self-adaptive  hp- 
FEM  are  comparable  (and  in  many  cases  more  accurate)  than  those  obtained  with  semi-analytical 
techniques  such  as  the  Mode  Matching  method,  for  problems  where  the  latest  methods  can  be 
applied. 


*The  authors  want  to  acknowledge  the  support  of  Ministerio  de  Education  y  Ciencia  of  Spain  under 
projects  TEC2010-18175/TCM  and  MTM2010-16511. 
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In  the  past  few  years,  there  have  been  increasing  research  activities  aiming  at  developing  a  robust  and 
efficient  numerical  solution  strategy  combating  multi-scale  electromagnetic  problems.  Such  multi-scale 
electromagnetic  problems  tax  heavily  on  numerical  methods  (finite  elements,  finite  difference,  integral 
equation  methods  etc.)  in  terms  of  desired  accuracy  and  stability  of  mathematical  formulations.  Besides, 
almost  various  physical  phenomena  are  usually  interacting  and  coupling  to  each  other.  For  instance,  the 
resistivity  of  most  conducting  metals  increases  linearly  with  the  increases  of  the  surrounding  temperature 
resulting  from  the  Joule  heating  by  electric  currents  flowing  through  the  conductors.  Therefore,  in  order 
to  accurately  characterize  the  performance  of  high  power  integrated  circuits  (ICs),  package  and  printed 
circuit  boards  (PCBs),  it  is  essential  to  account  for  both  electric  and  thermal  effects  and  the  intimate 
couplings  between  them. 


One  of  the  major  power  integrity  (PI)  issues  is  the  IR  drop,  caused  by  the  finite  resistivity  of  metals  and 
current  drawn  off  from  the  power/ground  planes.  When  the  resistance  of  power/ground  plane  increases 
with  the  shrinking  of  complex  geometries,  the  IR  drops  will  in  turn  increase  as  well.  To  further  confound 
the  difficult  issue,  the  rise  of  the  thermal  temperature  due  to  the  current  carrying  interconnects  also  has 
tremendous  impact  on  the  IC  performance  and  reliability.  Current  flow  in  a  VLSI  interconnect  can  cause 
a  power  dissipation,  which  is  referred  to  as  Joule  heating  or  self-heating.  Therefore,  due  to  the  linear 
temperature  dependency  of  metal  resistivity,  the  non-uniform  temperature  distribution  affects  the 
electrical  performance  of  power  delivery  network  (PDN)  and  substantially  increases  the  IR  drops  in  the 
power/ground  planes. 

In  this  paper,  we  present  non-conformal,  non-overlapping  domain  decomposition  methods  (DDMs)  for 
the  thermal-aware  DC  IR  drop  co-analysis  of  high-power  IC-Package-PCBs.  The  proposed  DDM  starts 
by  partitioning  the  composite  device  into  inhomogeneous  sub-regions  with  temperature  dependent 
material  properties.  Subsequently,  each  sub-domain  is  meshed  independently  according  to  its  own 
characteristic  features.  As  a  consequence,  the  troublesome  mesh  generation  task  for  complex  ICs  can  be 
greatly  subdued.  The  proposed  thermal-aware  DC  IR  drop  co-analysis  applies  the  non-conformal  DDM 
for  both  conduction  and  steady  state  heat  transfer  analyses  with  two-way  couplings  between  them. 
Numerical  examples,  including  an  IC  package  and  an  IC-Package-PCB,  demonstrate  the  flexibility  and 
potentials  of  the  proposed  thermal-aware  DC  IR-drop  co-analysis  using  non-conformal  DDMs. 
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In  this  paper,  we  propose  developing  a  hybrid  time  domain  discontinuous  Galerkin  vector  general¬ 
ized  finite  element  method  (TD-DG-VGFEM)  and  examine  its  accuracy,  efficiency,  and  stability  for 
transient  electromagnetic  (EM)  problems.  Both  VGFEM  and  DG  can  be  cast  under  the  umbrella 
of  partition  of  unity  methods.  In  VGFEM,  the  partition  of  unity  function  is  defined  on  overlapping 
domains,  whereas  in  DG,  the  overlap  is  zero.  As  a  consequence,  while  continuity  is  built  into  space 
of  functions  in  VGFEM,  it  is  not  in  DG.  Therefore,  additional  constraints  are  necessary  in  DG 
to  enforce  continuity.  Both  these  methods  have  been  used  to  provide  robust  solutions  to  electro¬ 
magnetic  problems  (O.  Tuncer,  B.  Shanker  and  L.C  Kempel,  IEEE  TAP  2010,  IEEE  APS  2011; 
J.  S.  Hesthaven  and  T.  Warburton,  J.  Comp.  Phys.,  2002;  I.  Perugia,  D.  Schotzau,  and  P.  Monk, 
Comput.  Methods  Appl.  Mech.  Engrg.  2002;  S.  Dosopoulos  and  Jin- Fa  Lee,  IEEE  TAP  2010). 
Both  methods  are  theoretically  capable  of  including  polynomials  and  non-polynomial  functions  in 
approximation  space.  While  non-polynomial  approximation  functions  have  been  used  in  VGFEM, 
they  have  not,  to  our  knowledge,  been  used  within  a  DG  framework. 

In  this  paper,  we  intend  to  investigate  a  hybridization  of  these  two  methods  to  solve  Maxwell 
equations  in  the  time  domain.  Specifically,  we  are  interested  in  examining  trade-offs  that  arise  out 
of  this  hybridization,  and  stability  and  convergence  of  the  proposed  approach.  To  this  end,  we  will 
present  the  mathematical  framework  necessary  for  this  integration  as  well  as  a  number  of  results 
that  examine  these  issues.  In  addition,  we  present  applications  to  a  range  of  practical  problems. 
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As  technological  devices  become  smaller,  faster,  and  more  complex,  computer  simulations  of  physical 
phenomena  become  an  increasingly  important  part  of  the  design  process.  The  multiscale  problem  caused 
by  a  disparate  mesh  usually  leads  to  a  large  unknown  count  and  a  slow  convergence  of  the  iterative  solver 
because  low-frequency  circuit  physics  and  high-frequency  wave  physics  coexist  in  one  problem. 
Applying  single  CEM  method  for  the  entire  problem  will  inevitably  cost  tremendously  more 
computational  resources  and  incapable  of  treating  every  part  of  the  problem  domain  efficiently.  For 
example,  it  is  most  efficient  to  solve  the  PEC  structures  with  boundary  integral  equation  (BIE)  accelerated 
with  powerful  fast  algorithms.  However,  the  complicated  dielectric  materials  are  more  efficient  to  be 
modeled  with  finite  element  method  (FEM).  For  the  eclectically  large  surrounding  regions,  high 
frequency  approximation  is  a  more  preferable  choice. 

Aimed  at  efficiently  calculating  the  radiation  and  scattering  from  arbitrary  composite  objects,  multiple 
targets  and  finite  periodic  structures,  a  multi-solver-based  generalized  impedance  boundary  condition 
(MS-GIBC)  is  proposed  and  finally  a  domain  decomposing  method  (DDM)  is  achieved.  According  to  the 
Huygens  equivalence  principle,  the  scattering  from  complex  structures  is  equivalent  to  the  radiation  of 
equivalent  sources  (electric  current  and  magnetic  current)  resided  on  their  boundaries.  A  generalized 
impedance  boundary  condition  (GIBC)  is  established  between  the  equivalent  magnetic  current  and 
electric  current.  Different  solvers,  such  as  BIE  and  FEM,  can  be  adopted  in  different  domain  to  construct 
the  GIBC.  Finally,  we  only  need  to  solve  the  boundary  integral  equations  concerning  equivalent  current  at 
each  boundary.  Moreover,  it  is  very  convenient  to  use  multilevel  fast  multipole  algorithm  (MLFMA), 
multifrontal  massively  parallel  sparse  direct  solver  (MUMPS)  and  various  preconditioners  to  accelerate 
the  solution  process. 

Due  to  the  GIBC,  only  equivalent  electric  current  Js  at  each  interface  is  introduced  and  a  sparse  matrix 
with  smaller  size  for  each  region  is  inverted  to  transfer  the  original  unknowns  to  the  unknowns  on  the 
equivalent  surface.  The  dimension  of  the  final  impedance  matrix  becomes  much  smaller  as  well. 
Furthermore,  the  non-conforming  grids  are  permitted  to  discretize  the  inner  fields  and  currents  on 
equivalent  surface.  Compared  with  recently  proposed  multi-solver  domain  decomposition  method  (MS- 
DDM),  the  tangential  continuities  of  E  and  H  fields  across  each  interface  are  imposed  directly  without 
extra  transmission  conditions.  Various  combinations  of  BIEs  in  our  method  right  correspond  to  different 
transmission  conditions  in  FEM-DDM. 
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In  this  paper,  the  dual-primal  finite  element  tearing  and  interconnecting  (FETI-DP)  method,  a  robust 
domain  decomposition  method,  is  adopted  together  with  the  LU  recombination  method  to  expand  the 
capability  and  improve  the  efficiency  of  3-D  finite-element  analysis  of  electromechanical  problems.  The 
FETI-DP  method  divides  an  original  large-scale  problem  into  smaller  subdomain  problems  and  then  deals 
with  the  subdomains  in  parallel  using  multiple  processors  to  significantly  reduce  the  total  computation 
time.  The  LU  recombination  method  is  adopted  to  deal  with  the  low-frequency  breakdown  problem, 
which  often  accompanies  the  finite-element  analysis  of  electromechanical  problems. 

A  general  electromechanical  device  consisting  of  linear  materials  can  be  described  by  the  magnetic 
diffusion  equation 

V  x  (vV  x  A)  +  jcocrA  =  J .  inflczT?3  (1) 

where  V  =  jU~l .  With  the  FETI-DP  method,  the  entire  computational  domain  is  first  partitioned  into  Ns 
subdomains.  At  each  subdomain  interface,  the  tangential  continuity  of  the  magnetic  field  intensity  and  the 
normal  continuity  of  the  magnetic  flux  density  should  be  satisfied,  which  can  be  written  as 

j«  x  ( h.  xAJ  =  hj  x  ( h  .  x  A ;. ) 

I «.  x (v. V x  A.)  =  -hj  x (VjV xA7)  =  A.  ^ 

The  second  boundary  condition  contains  an  unknown  vector  A ,  often  referred  to  as  the  “dual”  unknown, 
which  has  to  be  solved  to  recover  the  field  values  inside  each  subdomain.  In  order  to  enforce  the 
continuity  condition  along  corner  edges,  one  set  of  global  corner  unknowns,  Ac ,  often  referred  to  as  the 
“primal”  unknowns,  is  used  within  the  entire  computational  domain.  Application  of  the  standard  finite 
element  method  to  each  subdomain  yields  the  subdomain  system  which  involves  the  dual,  primal,  and 
remaining  unknowns.  By  eliminating  the  remaining  unknowns,  we  obtain  two  global  systems  of  equations 
involving  only  the  primal  and  dual  unknowns.  By  further  eliminating  the  primal  unknowns,  we  obtain  an 
interface  system  which  involves  only  the  dual  unknowns.  This  interface  system,  which  has  a  much 
smaller  number  of  unknowns  than  that  in  the  original  3-D  problem,  can  be  solved  by  a  Krylov  subspace 
method.  The  introduction  of  the  aforementioned  primal  unknowns  improves  the  convergence  rate  of  the 
Krylov  subspace  method  in  a  similar  manner  to  the  multigrid  method. 

When  dealing  with  electromechanical  problems,  the  finite-element  analysis  with  vector  basis  functions 
usually  encounters  the  so-called  low-frequency  breakdown  problem.  An  effective  way  to  solve  this 
problem  is  to  use  the  LU  recombination  algorithm  to  remove  the  dependant  equations  (H.  Ke,  T.  H. 
Hubing,  and  F.  Maradei,  “Using  the  LU  recombination  method  to  extend  the  application  of  circuit- 
oriented  finite  element  methods  to  arbitrarily  low  frequencies,”  IEEE  Trans.  Microw.  Theory  Tech.,  vol. 
58,  pp.  1189-1195,  2010).  With  the  LU  recombination  method,  the  system  matrix  is  first  factorized  to  be 
a  lower-  and  an  upper-triangular  matrices.  Since  a  null  space  exists  for  the  system  matrix,  some  diagonal 
values  of  the  upper-triangular  matrix  are  several  orders  smaller  than  others.  These  nearly-zero  diagonal 
values  actually  correspond  to  the  dependant  equations  of  the  system  matrix,  which  can  be  identified  easily. 
By  eliminating  the  dependant  equations  and  the  corresponding  unknowns,  the  system  matrix  is  then 
regularized.  Compared  to  the  tree-cotree  splitting  method,  this  method  deals  with  the  system  matrix 
directly,  and  therefore  can  be  implemented  more  easily  together  with  higher-order  basis  functions. 
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In  the  last  decade,  various  hybrid  electromagnetic-circuit  (EM-CKT)  simulators  that  simultaneously  solve 
Kirchoff  s  equations  for  voltages/currents  on  lumped-element  networks  and  Maxwell’s  equations  for 
electromagnetic  fields/currents  on  complex  volumes/surfaces  have  been  proposed.  On  the  one  hand,  such 
simulators  generally  use  the  modified  nodal  analysis  technique  to  formulate  circuit  equations  and  solve 
them  for  node  voltages  and  voltage-source  currents.  On  the  other  hand,  they  employ  a  variety  of 
techniques  to  formulate  partial  differential  or  integro-differential  electromagnetic  equations  and  solve 
them  for  fields  or  currents,  respectively,  and  can  be  categorized  as  differential-  or  integral-equation  based 
simulators.  These  hybrid  EM-CKT  simulators  can  also  be  categorized  based  on  their  temporal 
discretization  method  as  frequency-domain,  time-domain,  and  more  recently,  envelope-tracking  solvers 
(V.  Subramanian  and  A.  E.  Yilmaz,  “An  envelope-tracking  field-circuit  simulator  for  nonlinearly  loaded 
wire  antennas,”  in  Proc.  USNC/URSI  Radio  Science  Meeting,  July  2011);  they  represent  time  variations 
of  signals  (fields,  currents,  voltages)  of  interest  in  terms  of  sinusoidal  functions,  interpolator  sub-domain 
basis  functions,  and  Fourier  series  with  slowly  time-varying  coefficients,  respectively. 

This  article  focuses  on  integral-equation  based  time-domain  and  envelope-tracking  hybrid  EM-CKT 
simulators  and  their  coupling  methods.  The  principle  formulation  used  in  these  simulators  is  as  follows 
(A.  E.  Yilmaz,  J.  M.  Jin,  and  E.  Michielssen,  “A  parallel  FFT  accelerated  transient  field-circuit 
simulator,”  IEEE  Trans.  Microwave  Theory  Techniques,  vol.  53,  2005,  pp.  2851-2865):  In  the  CKT- 
solver  component,  first  the  vector  V(f),  which  represents  node  voltages  and  voltage-source  currents,  is 
discretized,  then  branch  constitutive  relations  at  voltage  source  branches  and  Kirchoff  s  current  law  at 
each  node  are  enforced,  i.e.,  Ik  =  0 ,  where  Ik  denotes  the  current  entering  to  the  node  from  the  kth 
branch  connected  to  it.  In  the  EM-solver  component,  first  the  current  densities  are  discretized,  then  the 
time-derivative  of  the  tangential  boundary  conditions  are  enforced,  e.g.,  dtE™n  (r,t)  =  -OtEZ  (J,  r,  t)  on 
a  perfectly  conducting  surface  S ,  where  J  denotes  the  surface  current  density  on  S  and  the  time- 
derivative  of  the  sum  of  components  of  the  incident  and  scattered  fields  tangential  to  S  are  equated  to 
zero.  The  two  components  are  interfaced  at  electrically  small  (lumped)  ports  using  controlled  sources: 
The  EM-solver  component  enforces  the  time-derivative  of  (line  integral  of)  the  total  electric  field  at  each 
port  p  be  equal  to  the  time-derivative  of  the  port  voltage  Vp(t),  which  is  discretized  by  the  CKT-solver 
component.  Similarly,  the  CKT-solver  component  enforces  the  currents  entering  to  the  nodes  of  each  port 
p  to  be  equal  to  that  of  the  current  leaving  the  port  JJ  J pds ,  which  is  discretized  by  the  EM-solver 
component.  While  this  J  /  dtE  —  V  /  I  formulation  is  suitable  for  coupling  established  EM  and  CKT 
solvers  without  major  changes  to  their  implementation,  it  requires  that  the  time-derivative  of  port  voltages 
be  computed  numerically,  which  can  cause  error  propagation  and  non-convergence  especially  when  the 
circuits  of  interest  contain  strong  nonlinearities.  In  this  article,  two  alternative  formulations  that  avoid  the 
numerical  time  derivative  in  the  coupling  are  proposed:  The  first  alternative  is  the  J  /  dtE  —  dt\  / 1 
formulation,  where  the  EM-solver  component  remains  the  same  as  above  but  the  CKT-solver  component 
is  reformulated  to  discretize  the  time-derivative  of  voltages.  The  second  alternative  is  the 
f  Jdt/E-V/  J  Idt  formulation,  where  the  EM-solver  component  is  reformulated  to  discretize 
integral  of  the  surface  current  density  and  to  enforce  the  tangential  boundary  conditions  (and  not  their 
time  derivatives)  and  the  CKT-solver  component  is  reformulated  to  enforce  the  time-integral  of  the 
Kirchoff  s  current  law  at  each  node.  Various  numerical  results  for  linearly  and  nonlinearly-loaded 
antennas  that  demonstrate  the  merits  of  these  alternative  formulations  will  be  presented  at  the  workshop. 
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Different  hybridizations  of  the  finite  element  method  (FEM)  and  the  method  of  moments  (MoM)  have 
been  extensively  utilized  for  an  exact  truncation  of  the  unbounded  spatial  domain  in  the  FEM  analysis  for 
quite  some  time.  As  one  of  the  principal  approaches  to  the  FEM-MoM  hybridization,  symmetric  coupling 
between  FEM  and  MoM  is,  essentially,  a  combination  of  the  appropriate  FEM  equations  (e.g.,  a  Galerkin- 
type  weak  form  of  the  curl-curl  electric  field  vector  wave  equation)  with  the  magnetic  field  integral 
equation  (MFIE)  and  the  electric  field  integral  equation  (EFIE)  that  gives  rise  to  a  symmetric  FEM-MoM 
matrix  system  of  equations. 

This  paper  presents  three  different  symmetric  and  nonsymmetric  FEM-MoM  techniques  using  higher 
order  basis  functions  and  generalized  curved  parametric  hexahedral  finite  and  quadrilateral  boundary 
elements,  with  both  FEM  and  MoM  leading  to  considerably  smaller  numbers  of  unknowns  when 
compared  to  similar  methods  with  low-order  basis  functions  and  simple  elements.  Curl-  and  divergence- 
conforming  hierarchical  polynomial  vector  bases  are  used  for  the  field  and  current  approximations  on 
FEM  hexahedra  and  MoM  quadrilaterals,  respectively,  enabling  different  coupling  methods  and  ways  of 
combining  appropriate  electromagnetic  quantities  while  satisfying  their  natural  interrelations. 

A  nonsymmetric  formulation  (M.  M.  Ilic,  M.  Djordjevic,  A.  Z.  Ilic,  and  B.  M.  Notaros,  IEEE  Trans. 
Antennas  Propagat.,  vol.  57,  May  2009)  places  equivalent  currents  at  the  boundary  surface  around  the 
FEM  domain  and  equates  tangential  electric  and  magnetic  field  components  on  each  side  of  the  boundary, 
while  not  enforcing  a  direct  connection  between  the  magnetic  current  and  the  electric  field  at  the  surface. 
It  uses  a  FEM  field  representation  in  the  MFIE  and  EFIE,  and  the  unknowns  are  surface  currents  at  the 
boundary,  which  are  later  used  to  calculate  the  fields  in  the  FEM  domain. 

In  a  higher  order  FEM-MoM  technique  developed  and  implemented  according  to  the  first  symmetric 
formulation  considered  (D.  J.  Hoppe,  L.  W.  Epp,  and  J.  Lee,  IEEE  Trans.  Antennas  Propagat.,  vol.  42, 
June  1994),  FEM,  EFIE,  and  MFIE  equations  are  used  independently  and  connected  in  a  final  system  of 
equations.  Again,  equivalent  surface  currents  are  placed  at  the  FEM  boundary,  giving  three  sets  of 
unknown  coefficients,  for  the  electric  field  and  the  magnetic  and  electric  currents,  respectively,  with  only 
an  implicit  connection  between  the  field  and  the  magnetic  current. 

In  another  higher  order  FEM-MoM  technique  developed  and  implemented  according  to  the  second 
symmetric  formulation  (M.  N.  Vouvakis,  S.  Lee,  K.  Zhao,  and  J.  Lee,  IEEE  Trans.  Antennas  Propagat., 
vol.  52,  Nov.  2004),  the  FEM  and  MFIE  are  combined  into  one  equation,  while  the  EFIE  is  independently 
added.  The  method  explicitly  forces  a  natural  relation  between  the  electric  field  and  the  magnetic  current, 
and  uses  the  same  coefficients  in  expansions  of  these  unknowns.  It  thus  computes  only  two  sets  of 
unknown  coefficients,  while  directly  solving  for  all  three  unknown  quantities.  In  support  of  this  option, 
the  method  implements  dual  sets  of  curl-  and  divergence-conforming  basis  functions  modeling  the  field 
and  current,  respectively. 

The  paper  evaluates,  discusses,  and  compares  theoretical  features  and  computational  properties  of  the 
three  methods,  all  implemented  using  truly  higher  order  geometrical  and  current/field  modeling. 
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Metamaterials  have  attracted  a  lot  of  attention  from  the  microwave  community  due  to  their 
unprecedented  electromagnetic  properties  and  novel  modes  (Special  Issue  on  Metamaterials,  IEEE 
Transaction  on  Antennas  and  Propagation,  Vol.  51,  No.  10,  2003,  pp.  2546-2750).  Many  new  antennas 
and  microwave  devices  have  been  developed,  harnessing  the  new  phenomena  realized  by  man-made 
metamaterials  such  as  negative  refraction.  The  latter  was  essential  in  proposing  sub-wavelength  focusing 
super  lenses  and  even  more  controversial  “cloaking”  applications.  We  recently  demonstrated  that  material 
anisotropy  provides  even  greater  mode  diversity.  Moreover,  inclusion  of  ferromagnetic  materials  leads  to 
more  degrees  of  design  freedom,  paving  way  for  unforeseen  microwave  devices  and  new  antenna  designs 
(J.L.  Volakis  et  al ,  “Narrowband  and  Wideband  Metamaterial  Antennas  Based  on  Degenerate  Band  Edge 
and  Magnetic  Photonic  Crystals,”  Proc.  IEEE,  Vol.  99,  No.  10,  2011,  pp.  1732-1745).  In  particular, 
biased  ferromagnetic  layers  used  in  a  simple  1  -dimensional  magnetic  photonic  crystal  (MPC)  can  support 
a  “frozen  mode”  with  zero  group  velocity.  Design  of  such  a  structure  using  the  finite  element  method  is 
perhaps  the  most  accurate  and  efficient  among  other  alternatives. 

Thanks  to  the  developments  in  the  past  decade,  the  FEM  is  quite  mature  in  modeling  infinite 
periodic  structures  as  well  as  antenna  arrays.  For  metamaterial  modeling,  periodic  boundary  conditions 
are  typically  utilized  to  truncate  the  FEM  domain  and  dispersion  characteristics  can  be  readily  computed. 
Geometrical  and  material  complexity  of  the  unit  cell  can  be  easily  handled  by  the  adaptability  of  the  FEM 
mesh.  However,  electromagnetic  parameters  of  the  constituent  materials  must  be  supplied  before  one  can 
form  the  FEM  eigensystem  and  its  subsequent  solution.  This  creates  an  issue  if  the  dielectric  materials 
have  frequency-dependent  permittivity  and/or  ferromagnetic  materials  have  permeabilities  as  a  function 
of  the  local  magnetic  field  intensity.  To  alleviate  this,  an  iterative  approach  can  be  employed  to  compute 
the  dispersion  characteristics  of  these  layouts.  Material  properties  of  the  substrate  and  the  computed 
dispersion  data  need  to  be  iteratively  updated  until  a  predefined  error  measure  is  satisfied. 

More  importantly,  experimental  verification  of  metamaterial  dispersion  characteristics  is  a 
challenge  due  to  finite  sizes  of  samples.  Evanescent/boundary  modes  usually  corrupt  measurement  data 
and  the  dispersion/group-delay  needs  to  be  indirectly  calculated.  For  this,  one  can  use  two  finite 
metamaterial  slabs  of  different  thicknesses  to  remove  the  evanescent  field  corruption  (N.  Apaydin  et  al ., 
’’Experimental  verification  of  frozen  mode  phenomenon  in  printed  magnetic  photonic  crystals,”  Proc.  5th 
European  Conf.  Antennas  and  Propagation,  2011,  pp.  2396  -2398).  We  should  remark  that  the  number  of 
unit-cells  in  the  finite  samples  should  be  chosen  to  be  reasonably  large,  such  that  the  coveted 
metamaterial  properties  are  prominent.  This  requirement  presents  an  inevitable  computational  challenge 
for  modeling  large,  finite-size  metamaterials  with  very  small  features  compared  to  the  wavelength.  We 
will  present  computational  results  to  illustrate  the  benefits  of  the  FEM  for  the  design  and  validation  of  this 
new  class  of  magnetic  metamaterials. 
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The  evaluation  of  fields  in  time  domain  for  large,  dense,  and/or  multiscale  periodic  structures  is 
a  topic  of  significant  interest  to  many  sectors  of  the  electromagnetics  community.  Characteristics 
such  as  rise-time,  ringing,  dispersive  effects,  time- varying  feeds  and  many  other  temporal  processes 
are  of  crucial  importance  to  accurate  designs  of  numerous  devices.  Additionally,  the  inclusion  of 
nonlinearities,  which  are  integral  to  many  electromagnetic  systems,  is  most  amenable  to  time  do¬ 
main  analysis.  Time  Domain  Integral  Equation  (TDIE)  solvers  are  a  promising  class  of  tools  for 
solving  periodic  systems  in  time  due  to  their  implicit  satisfaction  of  free  space  boundary  conditions 
without  requiring  domain  truncation  and  absorbing  boundary  conditions.  However,  the  evaluation 
of  temporal  and  spatial  convolutions  required  in  the  TDIE  kernel  is  a  computationally  intensive 
task,  scaling  as  0(7V27V2)  per  mode  when  fields  are  represented  in  terms  of  a  Floquet  expansion. 
Here,  Ns  is  the  number  of  spatial  degrees  of  freedom  and  Nt  is  the  number  of  time  steps.  Addition¬ 
ally,  the  Floquet  representation  of  the  time  domain  periodic  Green’s  function  consists  of  an  infinite 
summation  over  spectral  modes  that  must  be  truncated  in  some  way  to  be  practical. 


In  this  work  we  present  a  method  for  effecting  fast  convolutions  with  the  Time  Domain  Floquet 
Wave  Green’s  function  that  is  constructed  specifically  for  introduction  into  a  Marching  on  in  Time 
(MOT)  TDIE  framework.  To  accomplish  the  speedup,  we  utilize  an  alternate  representation  of  the 
Time  Domain  Floquet  Wave  Green’s  Function  that  explicitly  separates  space  and  time,  allowing 
independent  spatial  and  temporal  acceleration.  We  employ  the  method  of  Accelerated  Cartesian 
Expansions,  an  0{NS )  scheme  based  on  generalized  Taylor  expansions  to  accelerate  spatial  con¬ 
volutions,  and  a  block-FFT  scheme  to  accelerate  temporal  convolutions.  The  resulting  algorithm 
scales  as  0(NsNt  log2  Nt)  per  Floquet  mode.  In  addition  to  accelerating  the  convolutions,  we  also 
employ  a  strictly  band-limited  temporal  basis  function  which,  upon  convolution  with  the  time  do¬ 
main  Floquet  Wave  Green’s  function,  limits  the  number  of  contributing  Floquet  modes  to  0(1) 
(Chen  at  el.  Radio  Science  2005,  Vol  40).  The  form  in  which  the  temporal  convolutions  are  cast 
allows  reconstruction  to  arbitrary  precision  using  the  ACE  algorithm.  Thus,  the  entire  method 
scales  asymptotically  as  0(NsNt  log2  Nt)  and  can  represent  reconstructed  fields  to  a  high  degree  of 
accuracy.  In  this  work  we  provide  a  detailed  explanation  of  the  formulation,  bounds  on  error,  and 
results  validating  scaling  of  the  method. 
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Omnidirectional  light  propagation  in  two-dimensional  (2D)  photonic  crystals  (PCs)  is  investigated.  An 
efficient  approach  to  identifying  a  complete  photonic  band  gap  (PBG)  in  2D  PCs  has  been  developed.  The 
in-plane  band  structure  of  2D  photonic  crystals  is  calculated  by  an  adaptive  finite  element  method,  as 
implemented  in  FEMLAB.  With  adopting  the  suitable  boundary  conditions,  the  eigenvalues  can  be  easily 
and  rapidly  calculated  no  matter  how  complex  the  geometric  structures  are.  By  symmetry,  the 
omnidirectional  photon  density  of  states  (PDOS)  can  be  calculated  based  on  the  in-plane  dispersion 
relation  within  the  irreducible  Brillouin  zone.  The  PDOS  corresponding  to  both  the  radiation  and 
evanescent  waves  can  be  obtained  accurately  and  efficiently.  We  demonstrate  that  the  "complete  band 
gaps"  showed  by  some  previous  work  due  to  considering  only  the  radiation  modes  will  be  closed  by 
including  the  contributions  of  the  evanescent  modes.  Therefore,  the  complete  PBGs  do  not  exist  in  2D 
PCs  where  evanescent  waves  can  propagate.  These  results  are  of  general  importance  and  relevant  to  the 
spontaneous  emission  by  an  atom,  or  to  dipole  radiation  in  two-dimensional  periodic  structures. 
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Metamaterials  are  engineered  periodic  structures  whose  macroscopic  electromagnetic  response  is 
determined  by  the  subwavelength  features  of  their  constituent  unit  cells.  Accurate  characterization  of  this 
response  in  a  full-wave  simulation  environment  requires  dense  discretization  of  the  features  of  each  unit 
cell,  and  often  several  hundred  or  thousand  such  unit  cells  are  present  a  single  computational  domain. 
Whereas  the  availability  of  powerful,  low-cost  computers  has  enabled  the  solution  of  many  previously 
intractable  simulation  problems,  it  has  also  motivated  researchers  to  explore  increasingly  complex 
metamaterial  structures.  As  a  result,  full-wave  simulation  of  large  metamaterial  structures  is  still  plagued 
by  lengthy  simulation  times  and  excessive  memory  requirements.  Effective-medium  models  based  on 
simplified  analytical  descriptions  of  metamaterial  unit  cells  (e.g.  equivalent-circuit  models)  have  been 
successful  in  mitigating  these  issues;  however,  these  models  are  not  always  accurate,  and  this  is 
particularly  true  for  the  highly  complex  unit-cell  geometries  generated,  for  example,  from  fractals,  space¬ 
filling  curves,  and  genetic  algorithms.  Alternatively,  the  effective-medium  parameters  of  an  isotropic, 
homogeneous  metamaterial  can  be  extracted  by  inversion  of  the  simulated  scattering  parameters  of  a  thin 
slab  [D.  R.  Smith  et  al.  Physical  Review  5,  vol.  65,  no.  19,  p. 195104,  2002.].  This  procedure  is  suited  to 
full-wave  simulators  like  Ansys  HFSS,  which,  in  simple  rectangular  geometries,  can  apply  periodic 
boundary  conditions  to  limit  the  size  of  the  computational  domain.  However,  more  complex  geometrical 
arrangements  of  metamaterials,  such  as  cylindrical  or  spherical  arrangements,  are  not  so  easily  described 
by  such  simulators,  particularly  because  these  arrangements  may  result  in  anisotropic  and  inhomogeneous 
properties.  The  computational  domains  for  these  structures  cannot  easily  be  reduced. 

To  address  this  issue,  we  present  a  hybrid  numerical/semi-analytical  approach  to  the  modeling  of 
generally  anisotropic  and  inhomogeneous  metamaterials.  Inhomogeneity  is  modeled  in  a  discrete, 
multilayer  fashion.  Anisotropy  is  addressed  by  decomposing  propagation  into  orthogonal  directions,  and 
an  equivalent  full-wave  rectangular  metamaterial  geometry  is  determined  for  each  direction.  Each  model 
may  then  be  reduced  using  periodic  boundary  conditions  and  effective-medium  parameters  corresponding 
to  each  direction  may  then  be  extracted  as  for  homogeneous,  isotropic  metamaterial  slabs.  An  analytical 
formulation  utilizes  these  extracted  parameters  in  an  effective  medium  model  to  obtain  the  near-  and  far- 
field  of  the  structure.  We  focus  on  the  case  of  a  cylindrical  negative-refractive-index  transmission-line 
(NRI-TL)  metamaterial  structure  [Pollock  et  al.  Proc.  of  Metamaterials  201 7,  30-32  (2011)]  created  by 
radially  arranging  several  planar  NRI-TL  metamaterial  layers,  which  exhibits  uniaxial  anisotropy  and 
radial  inhomogeneity.  We  apply  the  extracted  data  to  our  analytical  formulation  and  demonstrate  that  the 
near-fields,  radiation  patterns,  and  radiated  power  produced  by  a  source  placed  next  to  the  cylindrical 
NRI-TL  metamaterial  are  in  excellent  agreement  with  those  acquired  from  full-wave  FEM  simulations  of 
the  entire  structure  using  Ansys  HFSS.  The  full-wave  simulation  region  consisted  of  a  cylinder  of 
diameter  between  0.66  and  0.73  free-space  wavelengths  over  the  simulated  frequency  range,  and  a  height 
of  between  0.055  and  0.062  free-space  wavelengths.  A  sufficiently  converged  result  required  a  minimum 
of  990,000  tetrahehdra,  over  1 1GB  of  memory,  and  1.5  hours  of  simulation  time  per  frequency  point  on  a 
3.4GHz  quad-core  desktop  computer.  In  our  analytical  method,  the  solution  is  based  on  a  cylindrical- 
wave  expansion  of  the  fields,  which  causes  the  computation  time  to  depend  on  both  the  required  spatial 
resolution  and  truncation  limit  on  the  number  of  azimuthal  modes  contributing  to  the  solution.  To  obtain  a 
commensurate  degree  of  resolution  and  convergence,  the  fields  can  be  obtained  for  a  single  frequency 
point  in  just  over  5  minutes.  Even  when  taking  into  account  the  time  required  to  extract  the  anisotropic 
parameters  using  full-wave  simulations,  we  observe  an  order  of  magnitude  reduction  in  total  simulation 
time. 
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The  dual-primal  finite  element  tearing  and  interconnecting  (FETI-DP)  method  is  one  of  the  most 
advanced  non-overlapping  domain  decomposition  methods  (DDMs).  It  has  been  successfully  used  in 
large-scale  electromagnetic  analysis  such  as  phased-array  antenna  and  photonic  crystal  (PhC)  simulations. 
The  earlier  two  conformal  versions  of  this  method  in  computational  electromagnetics  community  are 
referred  to  as  FETI-DPEM1  (Y.  J.  Li  and  J.  M.  Jin,  “A  vector  dual-primal  finite  element  tearing  and 
interconnecting  method  for  solving  3-D  large-scale  electromagnetic  problems,”  IEEE  Trans.  Antennas 
Propagat .,  vol.  54,  pp.  3000-3009,  2006)  and  FETI-DPEM2  (Y.  J.  Li  and  J.  M.  Jin,  “A  new  dual-primal 
domain  decomposition  approach  for  finite  element  simulation  of  3D  large-scale  electromagnetic 
problems,”  IEEE  Trans.  Antennas  Propagat .,  vol.  55,  pp.  2803-2810,  2007),  respectively.  These  two 
schemes  employ  the  Dirichlet-Newmann  transmission  condition  and  the  Robin-type  transmission 
condition  at  the  subdomain  interfaces,  respectively.  Eigenspectrum  analysis  and  Numerical  simulation 
show  the  FETI-DPEM2  outperforms  FETI-DPEM1  in  terms  of  convergence  rate  of  the  iterative  solution 
of  the  global  interface  problem,  even  without  preconditioning.  At  the  same  time,  a  nonconformal 
counterpart  of  the  conformal  FETI-DP  method,  which  is  called  FETI-like  method,  was  proposed  (K.  Zhao 
et  al,  “A  domain  decomposition  method  with  nonconformal  meshes  for  finite  periodic  and  semi-periodic 
structures,”  IEEE  Trans.  Antennas  Propagat .,  vol.  55,  pp.  2559-2570,  2007).  It  reformulates  the  Robin- 
type  transmission  condition  and  has  a  similar  convergence  behavior  as  that  of  FETI-DPEM2. 

In  order  to  further  enhance  the  capability  of  the  FETI-DP  method  for  modeling  large-scale  and  complex 
antenna  structures,  two  nonconforming  FETI-DP  methods  are  proposed  in  our  recent  work.  They  are 
referred  to  as  the  nonconformal  FETI-DP  method  and  the  cement  FETI-DP  method,  respectively.  Both  of 
these  two  methods  implement  a  Robin-type  transmission  condition  at  the  subdomain  interfaces  to 
preserve  the  fast  convergence  when  solving  the  global  interface  problem  at  high  frequencies  iteratively. 
Similar  to  the  conformal  FETI-DP  method,  a  global  coarse  problem  related  to  the  unknowns  associated 
with  the  subdomain  comer  edges  is  designed  to  propagate  the  residual  error  to  the  whole  computational 
domain  at  each  iteration.  The  nonconformal  FETI-DP  method  extends  the  conformal  FETI-DP  method  to 
deal  with  nonconformal  interface  and  corner  meshes,  and  the  cement  FETI-DP  method  formulates  the 
FETI-like  method  to  combine  interface  dual  unknowns  with  global  primal  unknowns  defined  on  corners. 

The  FETI  method  is  particularly  suitable  for  simulating  large  finite  arrays  because  the  geometrical 
repetitions  of  an  array  structure  can  be  exploited  to  further  enhance  the  computational  performance.  In 
some  engineering  applications,  we  may  encounter  stmctures  with  repetitions  in  all  three  dimensions.  For 
example,  negative  refraction  by  using  materials  with  nonnegative  permittivity  and  permeability  is  a  hot 
research  topic  in  the  optics  society.  In  the  past  few  decades,  experiments  and  simulations  have 
demonstrated  negative  refraction  and  negative  refraction  imaging  by  two-dimensional  (2D)  PhC  flat 
lenses.  However,  these  2D  flat  lenses  have  shown  lensing  ability  only  along  one  direction;  consequently, 
the  usage  of  2D  flat  lenses  is  limited  in  real-life  applications.  For  this  reason,  “flat  spherical”  lenses, 
which  work  by  full  3D  negative  refraction  (F3DNR)  and  provide  super  resolution  imaging,  are  especially 
desired.  In  a  recently  published  paper,  a  so-called  3D  photonic  crystal  flat  lens  was  fabricated  and  the 
F3DNR  phenomenon  was  experimentally  observed.  This  paper  will  focus  on  simulating  the  F3DNR 
phenomenon  and  show  the  effectiveness  and  efficiency  of  the  nonconformal  and  cement  FETI-DP 
methods  in  modeling  3D  photonic  crystal  array  structures. 
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It  is  a  well  established  that  the  hybrid  finite  element  -  boundary  integral  -  multilevel  fast  multipole 
algorithm  (FE-BI-MLFMA)  is  a  general,  accurate,  and  efficient  computation  technique  for  open  region 
problems  such  as  scattering  /  radiation  problems.  The  conventional  algorithm  of  FE-BI-MLFMA  is 
directly  to  apply  the  iterative  solvers  to  the  final  FE-BI  matrix  equation  (X.  Q.  Sheng  et  al.  “On  the 
formulation  of  hybrid  finite-element  and  boundary-integral  method  for  3D  scattering,”  IEEE  Trans. 
Antennas  Propagat ,  Vol.  46,  No.  3,  1998,  pp.  303-311).  This  conventional  algorithm  converges  slowly, 
sometimes  even  failing  to  converge  due  to  ill-conditioned  nature  of  the  associated  FE-BI  matrix  equation. 
Many  fast  algorithms  of  FE-BI-MLFMA  have  been  developed  over  the  past  years  to  address  this 
shortcoming  (J.  Liu  et  al.,  “A  highly  effective  preconditioner  for  solving  the  finite  element-boundary 
integral  matrix  equation  for  3-D  scattering,”  IEEE  Trans.  Antennas  Propagat.  Vol.  50,  No.  9,  2002,  pp. 
1212-1221).  The  key  to  designing  fast  algorithms  of  FE-BI-MLFMA  is  the  efficient  construction  of  a 
suitable  preconditioner  for  the  FE-BI  matrix  equation.  Thereby,  the  FE-BI  matrix  is  first  decomposed  into 
parts  taking  advantage  of  inherent  favourable  properties  of  FEM  and  BI.  Different  methods  have  been 
employed  to  construct  preconditioners  for  the  FEM-  and  Bi-parts,  respectively.  In  particular  the  idea  of 
the  finite  element  tearing  and  interconnecting  method  (Y.  J.  Li  et  al.,  “A  new  dual-primal  domain 
decomposition  approach  for  finite  element  simulation  of  3-D  large-scale  electromagnetic  problems”, 
IEEE  Trans.  Antennas  Propag.,  vol.  55,  No.  10,  2007,  pp.  2803-2810)  is  employed  to  construct  the 
preconditioner  for  the  FEM  part  of  the  FE-BI  matrix  equation.  The  sparse  approximation  inverse  method 
is  directly  utilized  to  construct  the  preconditioner  for  the  BI  part  of  the  FE-BI  matrix  equation.  Afterwards, 
a  solution  algorithm  of  FE-BI-MLFMA  based  on  the  above  two  preconditioners  is  designed.  The 
numerical  performance  of  the  proposed  domain-decomposition-based  preconditoner  of  FE-BI-MLFMA 
is  investigated  in  detail,  including  the  scalability  of  this  algorithm. 
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In  computational  electromagnetics  and  other  application  fields,  multifrontal  method  is  considered  as  one 
of  the  most  efficient  direct  solvers  for  large  sparse  equation  system.  But  like  other  direct  solvers,  the 
multifrontal  method  requires  complete  construction  of  the  global  stiffness  matrix,  which  may  be  a 
significant  demerit  for  large-scale  applications.  Because  in  the  original  frontal  method,  the  global  stiffness 
matrix  is  never  constructed,  each  stage  of  assembling  and  solving  stiffness  matrix  is  alternately  executed, 
thus  fully  assembled  unknowns  are  eliminated  immediately  from  the  frontal  matrix.  This  process  makes 
the  frontal  method  based  programs  particularly  memory  efficient  for  large-scale  applications. 

To  overcome  aforementioned  demerit,  the  idea  of  partitioning  the  whole  computational  mesh  into  a 
number  of  subdomains  and  then  implementing  the  frontal  scheme  separately  in  each  subdomain  has  been 
appeared.  This  method  was  originally  termed  as  the  multiple  front  method.  Fully  implemented  multilevel 
domain-wise  multifrontal  approach  is  proposed  by  Kim:  (J.  H.  Kim  and  S.  J.  Kim,  “A  Multifrontal  Solver 
Combined  with  Graph  Partitioners,”  AIAA  Journal,  Vol.  38,  1999,  pp.  964-970),  where  a  bisection 
algorithm  is  first  applied  to  recursively  partition  the  sets  of  elements  into  two  subdomains,  and  then  two 
neighboring  domains  are  merged  in  the  reverse  order  until  all  elements  are  assembled  on  the  coarsest 
mesh.  In  these  approaches,  domain  partitioning  and  numbering  are  performed  instead  of  matrix  ordering, 
the  global  assembly  of  the  stiffness  matrix  is  not  required,  and  the  out-of-core  technique  can  be  used 
easily  and  seamlessly  as  the  frontal  method  does. 

Although  various  implementations  of  the  domain-wise  multifrontal  approach  have  been  studied  over  the 
years,  the  theoretical  knowledge  of  how  they  perform  lags  way  behind.  In  this  work,  the  computational 
complexity  of  a  P-way  domain  wise  multilevel  multifrontal  algorithm,  pertaining  to  the  finite  element 
method,  is  analyzed  on  regular  2D  and  3D  meshes  in  various  aspects  and  with  different  parameters,  such 
as  the  number  of  fronts,  the  number  of  partition  levels,  the  ordering  of  fronts,  and  the  type  of  boundary 
condition. 

The  operation  counts  of  both  local  and  multilevel  global  condensation  phase  are  accurately  estimated 
along  with  the  relevant  parameters,  which  make  the  theory  more  accurate  and  more  generous  than 
previous  estimations,  and  more  over,  enable  us  to  illustrate  how  the  number  of  fronts,  the  number  of 
partition  levels,  and  the  performance  of  the  method  are  related,  and  how  the  multilevel  multifrontal 
scheme  reduces  the  total  numerical  operations  intrinsically. 

The  operation  count  and  execution  time  of  the  method  are  measured  experimentally  on  a  sequence  of  2D 
and  3D  examples.  The  theoretical  estimations  agree  well  with  experimental  results,  as  well  as  the  optimal 
number  of  fronts.  The  results  also  show  that  due  to  its  intrinsical  reduction  of  the  operation  counts,  the 
domain  wise  multilevel  multifrontal  method  achieves  a  huge  amount  of  CPU  time  saving,  and  thereby 
offers,  without  approximation,  an  easy  way  to  substantially  upgrade  the  frontal  method  based  programs. 
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When  high  power  densities  are  involved  within  passive  microwave  and  millimeter-wave  devices,  materials 
may  present  a  significant  nonlinear  behavior.  A  frequency  domain  analysis  is  sometimes  preferable  to  a  time- 
domain  analysis,  especially  when  steady-state  fields  are  aimed  at  for  a  finite  set  of  frequencies.  Frequency 
domain  finite  element  (FE)  method  provides  an  invaluable  analysis  tool  for  linear  devices,  but  it  fails  when 
material  properties  present  a  dependence  on  the  fields  strength.  To  take  into  account  of  harmonic  distortion 
that  occurs  in  such  a  case,  Harmonic  Balance  (HB)  techniques  may  be  used. 

The  proposed  method  extends  the  previous  work  (G.  Guarnieri,  G.  Pelosi,  L.  Rossi,  and  S.  Selleri,  “A 
Domain  Decomposition  Technique  for  Efficient  Iterative  Solution  of  Nonlinear  Electromagnetic  Problems,” 
IEEE  Trans.  Antennas  Propag.,  vol.  58,  no.  12,  pp.  4090-4095,  Dec.  2010)  to  include  higher  order 
harmonics.  As  a  result,  a  global  HB-FE  system  is  assembled  and  solved  iteratively,  converging  to  the 
solution  within  a  prescribed  tolerance.  Depending  on  the  nonlinearity  degree  of  the  model,  either  Picard  or 
relaxed  iterations  are  used. 

One  or  two  fundamental  frequency  excitations  have  been  investigated,  the  latter  allowing  to  accurately 
compute  third  order  intermodulation  products.  Coupling  between  harmonics  is  performed  by  testing  the 
product  of  the  nonlinear  material  property  with  the  multiharmonic  field  amplitude  (F.  Bachinger,  U.  Langer, 
and  J.  Schoberl,  “Numerical  analysis  of  nonlinear  multiharmonic  eddy  current  problems,”  Numer.  Math., 
vol.  100,  no.  4,  pp.  593-616,  Jun.  2005).  The  material  property  can  actually  be  expressed  in  terms  of  a 
sum  of  harmonics,  being  dependent  on  the  multiharmonic  field.  The  procedure  of  testing  can  be  performed 
analytically  for  the  single  frequency  excitation  case,  providing  algebraic  relations  between  first  and  higher 
order  harmonics  (S.  Yamada,  K.  Bessho,  and  J.  Lu, “Harmonic  balance  finite  element  method  applied  to 
nonlinear  AC  magnetic  analysis,”  IEEE  Trans.  Magn.,  vol.  25,  no.  4,  pp.  2971-2973,  Jul.  1989).  In  case  of 
multiple  excitation  frequencies  impinging  on  the  device,  an  analytic  testing  can  become  impractical.  Hence, 
coupling  coefficients  have  been  computed  by  the  use  of  the  fast  Fourier  transform  algorithm. 

A  substructuring  Domain  Decomposition  method  has  been  applied  to  isolate  unknowns  contributing  to  non¬ 
linear  elements  from  those  remaining  linear,  in  order  to  limit  the  nonlinear  solving  loop  to  a  smaller  domain. 
Numerical  simulations  have  been  carried  out  to  validate  the  technique  and  perform  timings  comparisons 
with  the  full  domain  HB-FE  analysis. 
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The  dual-primal  finite  element  tearing  and  interconnecting  (FETI-DP)  method,  a  robust  domain 
decomposition  method,  has  been  successfully  applied  to  solve  large-scale  electromagnetic  problems  (Y. 
Li  and  J.  M.  Jin,  “A  vector  dual-primal  finite  element  tearing  and  interconnecting  method  for  solving  3-D 
large-scale  electromagnetic  problems,”  IEEE  Trans.  Antennas  Progag .,  vol.  54,  no.  10,  pp.  3000-3009, 
2006).  The  FETI-DP  method  divides  an  original  large-scale  problem  into  smaller  subdomain  problems 
and  then  deals  with  the  subdomains  in  parallel  using  multiple  processors  to  significantly  reduce  the  total 
computation  time.  However,  when  the  operating  frequency  is  very  low,  the  finite-element  analysis  with 
vector  basis  functions  results  in  an  ill-conditioned  system  matrix,  which  causes  the  so-called  low- 
frequency  breakdown  problem.  This  problem  may  also  occur  when  very  small  elements  are  used  in  the 
finite-element  mesh.  In  this  paper,  the  tree-cotree  splitting  (TCS)  method  is  adopted  to  eliminate  the  low- 
frequency  breakdown  problem  within  the  FETI-DP  framework. 

The  proposed  algorithm  is  used  to  solve  the  vector  Helmholtz  equation  for  the  electric  field  given  by 

V  x  (//-1 V  xE)  -  co2  =  -jcoii  inOcT?3  (1) 

where  £  and  /a  are  the  permittivity  and  permeability,  respectively.  The  FETI-DP  method  generally 
consists  of  three  steps.  First,  the  computational  domain  is  divided  into  Ns  nonoverlapping  subdomains, 
and  an  incomplete  solution  of  each  subdomain  is  obtained  using  a  direct  solver.  Next,  Dirichlet  and 
Neumann  boundary  conditions  are  enforced  at  the  subdomain  interfaces  using  Lagrange  multipliers, 
which  results  in  a  global  interface  and  a  global  comer  system.  These  two  global  systems  are  then  solved 
using  the  Krylov  subspace  method,  whose  solution  is  then  used  to  recover  the  field  values  inside  each 
subdomain.  The  constmction  of  the  global  corner  system  improves  the  convergence  rate  of  the  Krylov 
subspace  method  in  a  similar  manner  to  the  multigrid  method.  In  all  three  steps  of  the  FETI-DP  method, 
computations  can  be  accelerated  using  parallel  computing  schemes  to  reduce  the  computation  time. 

When  the  operating  frequency  is  very  low,  the  displacement  current  term  becomes  much  smaller 
compared  to  the  curl-curl  term  in  (1).  Unfortunately,  the  finite  element  discretization  of  the  curl-curl 
operator  with  edge  elements  result  in  an  ill-conditioned  matrix,  which  causes  the  low-frequency 
breakdown  problem.  An  effective  approach  to  solving  this  problem  is  to  use  the  TCS  method  to  remove 
the  dependent  equations  (J.  B.  Manges  and  Z.  J.  Cendes,  66 A  generalized  tree-cotree  gauge  for  magnetic 
field  computation,”  IEEE  Trans.  Magn .,  vol.  31,  pp.  1342-1347,  1995).  To  apply  the  TCS  method  in  the 
FETI-DP  framework,  three  conditions  must  be  satisfied  simultaneously:  (1)  A  tree  structure  is  constructed 
for  the  global  corner  system  consisting  of  all  the  corner  edges;  (2)  A  tree  structure  is  constructed  for  each 
subdomain  consisting  of  local  volume  and  interface  edges;  (3)  The  tree  and  cotree  edges  are  defined 
consistently  on  the  subdomain  interfaces  so  that  the  Dirichlet  and  Neumann  boundary  conditions  can  be 
enforced.  In  order  to  satisfy  all  three  requirements,  we  apply  the  TCS  algorithm  in  a  nested  manner.  First, 
we  construct  a  tree  structure  for  the  corner  system  to  connect  all  the  corner  nodes,  which  together  with  the 
original  reference  nodes  are  considered  as  the  new  root.  Next,  we  grow  the  tree  from  the  new  root  to 
include  all  the  nodes  on  the  subdomain  interfaces,  so  that  a  consistent  tree  structure  can  be  built  on  the 
interfaces.  Once  all  the  corner  and  interface  nodes  are  connected,  we  continue  to  construct  the  tree 
structure  inside  each  subdomain  until  a  complete  tree  is  formed. 
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For  solving  scattering  problems  of  polygonal  prisms,  the  normal  component  of  the  magnetic  field  may 
have  a  discontinuity  across  a  material  interface,  and  its  direction  changes  infinitely  rapidly  at  a  sharp  edge. 
The  node-based  finite  element  method  cannot  model  these  effects  properly  because  its  basis  functions  are 
continuous  in  all  components  from  one  element  to  other.  Fortunately,  edge  elements  are  more  flexible. 
Generally,  vector  fields  which  are  represented  by  edge  elements  have  normal  discontinuities,  and  are  only 
tangentially  continuous.  In  addition,  the  degrees  of  freedom  of  edge  elements  are  associated  with  the 
edges,  not  the  nodes,  so  they  allow  the  vector  field  to  change  direction  abruptly  at  the  sharp  edge. 
Therefore,  edge  elements  are  very  suitable  for  solving  scattering  from  polygonal  prisms. 

When  the  scattering  problem  is  discretized  by  edge  elements,  the  problem  is  transformed  into  a  linear 
system  of  algebraic  equations.  Generally,  electrically  large  objects  involve  a  large  number  of  unknowns, 
because  there  must  be  enough  edge  elements  per  wavelength  to  represent  the  field  variation  accurately. 
The  corresponding  linear  system  of  algebraic  equations  is  very  costly  to  solve.  For  such  systems,  the  idea 
of  splitting  the  domain  is  attractive.  The  whole  domain  is  divided  into  many  smaller  sub-domains,  which 
are  solved  individually.  The  size  of  each  sub-domain  can  be  adjusted  according  to  the  computer  resources 
available.  The  Decomposition  Projective  Method  (DPM)  can  meet  these  requirements  well  and  is 
therefore  used  here.  Not  only  is  the  order  of  the  system  of  equations  on  each  sub-domain  very  small,  but 
also  the  convergence  rates  of  DPM  is  geometric. 

With  DPM,  the  equations  on  sub-domains  are  solved  repeatedly  because  DPM  is  an  iterative  method. 
When  most  sub-domains  have  the  same  structure,  their  coefficient  matrix  of  the  linear  systems  of 
algebraic  equations  are  also  the  same.  If  the  basic  set  of  solutions  of  these  equations  is  solved,  the  solution 
can  be  got  by  combination  of  the  basic  set  of  solutions,  which  saves  much  time  on  computation.  Partial 
basic  solution  vectors  (PBSV)  method  can  complete  it. 

As  a  numerical  example,  edge  elements  and  DPM  with 
PBSV  are  used  to  analyze  scattering  from  a  square, 
perfectly-conducting,  infinitely-long,  cylinder  of  side 
10M).  Fig.  1  is  the  surface  current  density  on  the  surface 
of  the  cylinder  with  an  incident  TM  wave.  The  mesh 
nodes  are  almost  the  same  for  edge  elements  and  MoM. 

There  are  about  52,600  unknowns  in  the  edge  element 
method.  DPM  without  PBSV  takes  181s  to  compute,  and 
there  are  830  iterations.  Using  DPM  with  PBSV,  it  takes 
51s  and  also  830  iterations.  Clearly,  DPM  with  PBSV 
costs  less  time  than  DPM  alone.  Meanwhile,  Fig.  1  also 
shows  that  the  numerical  approach  of  edge  elements  is 
better  than  that  of  MoM  especially  along  side  CD. 


Fig.l.  Surface  current  density  on  the  cylinder. 
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Previously,  the  spectral  element  method  was  implemented  to  calculate  the  transmission  rate  of  an 
infinite  periodic  photonic  crystal  slab  under  plane  wave  incidence.  The  electric  field  in  the  com¬ 
putational  domain  is  interpolated  by  the  spectral  element  basis  functions.  The  periodic  boundary 
conditions  are  used  for  the  vertical  boundary,  and  the  method  of  moments  are  used  as  the  boundary 
conditions  for  the  top  and  bottom  horizontal  boundary. 

In  this  work,  the  implementation  was  extended  to  calculate  the  eigenstate  of  the  photonic  crystal 
slab.  The  method  of  moments  gives  integral  equation  of  electromagnetic  field  base  on  the  Green’s 
function,  which  is  a  nonlinear  function  of  the  frequency.  Because  the  frequency  is  the  eigenvalue,  the 
weak  form  of  the  spectral  element  method  with  the  method  of  moments  gives  a  nonlinear  eigenvalue 
problem  for  the  eigenstate  of  the  photonic  crystal  slab.  This  nonlinear  eigenvalue  problem  was 
solved  by  iterative  solver.  The  band  structure  of  the  air  bridge  photonic  crystal  can  be  obtained. 

In  addition,  the  second  order  nonlinear  optical  effect  is  also  modeled  in  this  numerical  method. 
With  the  given  nonlinear  susceptibility  of  the  nonlinear  optical  material,  the  second  order  nonlinear 
optical  effect  can  be  directly  simulated  by  the  Maxwell’s  wave  equation  of  the  electromagnetic  field. 
The  only  approximation  in  the  modeling  is  the  numerical  discretization  by  the  spectral  element 
method.  The  nonlinear  polarization  field  can  be  calculated  from  the  interpolation  of  the  electric 
field.  This  polarization  field  is  tested  as  the  source  field  in  the  spectral  element  method,  which  gives 
the  second  harmonic  generation  field.  The  second  harmonic  generation  rate  versus  the  wavelength 
and  incident  angle  were  calculated. 

The  band  structure  of  the  air  bridge  photonic  crystal  slab  is  compare  with  the  second  harmonic 
generation  rate.  Numerical  results  show  that  the  resonance  between  the  band  structure  and  the 
second  harmonic  frequency  could  enhance  the  nonlinear  effect.  The  additional  conditions  for  the 
enhancement  is  the  matching  of  the  symmetry  property  between  the  eigenstate  and  the  generated 
second  harmonic  polarization  field.  Comparing  to  homogeneous  dielectric  slab,  the  air  bridge 
photonic  crystal  can  enhance  the  second  harmonic  generation  rate  for  four  magnitude  at  certain 
wavelength. 
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Recently,  the  domain  decomposition  methods  (DDMs)  have  been  demonstrated  effective  in  solv¬ 
ing  many  mission-critical  electromagnetic  applications.  One  category  of  the  DDMs  is  the  non¬ 
overlapping  DDMs.  The  basic  idea  is  to  decompose  the  original  entire  problem  domain  into 
several  non-overlapping  sub-domains.  The  continuities  of  electromagnetic  fields  at  the  interfaces 
between  adjacent  sub-domains  are  enforced  through  some  boundary  conditions.  Moreover,  the  non- 
conformal  DDMs  which  relieve  the  restriction  of  mesh  conformity  have  been  proposed  and  shown 
to  be  accurate  and  efficient. 

We  present  a  non-overlapping  and  non-conformal  domain  decomposition  method  (DDM)  for  solving 
the  time-harmonic  Maxwell  equations  in  M3.  It  has  been  recognized  in  the  past  that  the  conver¬ 
gence  behavior  of  non-overlapping  DDMs  depends  heavily  upon  the  transmission  conditions  used 
to  enforce  the  continuity  of  tangential  fields  on  the  interface  between  sub-domains.  Recently,  a 
true  second  order  transmission  condition,  which  involves  two  second-order  transverse  derivatives,  is 
proposed  and  allows  convergence  of  both  propagating  and  evanescent  electromagnetic  waves  across 
domain  interfaces.  Here,  we  introduce  a  new  optimal  second  order  transmission  condition  to  further 
improve  the  convergence  in  the  DDM  iterations.  Furthermore,  through  an  analysis  of  the  DDM 
with  the  SOTC,  we  show  that  there  still  exists  a  weakly  convergent  region  where  the  convergence 
in  the  DDM  can  still  be  unbearably  slow  for  electrically  large  problems.  It  is  found  that  the  weakly 
convergent  region  is  centered  at  the  cutoff  modes,  or  electromagnetic  waves  propagate  in  parallel  to 
the  domain  interfaces.  Subsequently,  a  global  plane  wave  deflation  technique  is  utilized  to  derive  an 
effective  global-coarse-grid  preconditioner  to  promote  fast  convergence  of  the  cutoff  or  near  cutoff 
modes  in  the  vicinity  of  domain  interfaces. 
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Proliferation  of  CPU  cores  in  modern  computing  has  renewed  the  interest 
in  computational  strategies  and  algorithms  that  leverage  parallel  processing. 
In  the  area  of  scientific  computation,  such  as  computational  electromagnetics 
(CEM),  the  domain  decomposition  (DD)  computational  paradigm  offers  a  vi¬ 
able  alternative  to  traditional  finite  element  or  boundary  element  fast  solvers, 
e.g.  multigrid  or  fast  multipole.  DD  methods  first  divide  the  discrete  com¬ 
putational  model  into  portions  (domains)  that  are  individually  solved  in  each 
core.  During  the  second  and  most  important  step,  the  solution  of  the  overall 
model  is  recovered  by  locally  communicating  portions  from  each  domain  so¬ 
lution  across  processors,  often  in  an  iterative  fashion.  Although  this  strategy 
offers  distinct  advantages,  such  as  straightforward  programing  and  balanced 
usage  of  processing  power,  when  difficult  CEM  models  are  encountered,  no 
guarantees  on  speed  or  success  of  the  iterative  convergence  exist.  The  iter¬ 
ative  convergence  of  most  popular  DD  methods,  such  as  the  finite  element 
tearing  and  interconnecting  (FETI)  method,  depends  strongly  on  discretiza¬ 
tion  size,  frequency,  field  regularity  and  domain  shape,  hampering  scalability 
and  robustness. 

Over  the  past  years  a  large  body  of  literature  has  been  published  in  applied 
mathematics  and  CEM  communities  trying  to  overcome  the  various  robustness 
and  scalability  problems  using  higher  order  transmission  conditions  (HOTCs) 
(Rawat  et  al.  SIAM  J.  Sci.  Comput.  Vol.  32,  p.  3584,  2010),  combination  of 
primal  FEMs  with  FETI  (FETI-DP)  (Toselli,  IMA  J.  Numer.  Anal.,  Vol.  26, 
p.  96,  2006),  (Y.  Li  et  ah,  IEEE  Trans.  Antennas  Prop.,  Vol.  55,  p.  2803, 
2007),  or  corner  penalty  methods  (Z.  Peng  et  al.  IEEE  Trans.  Antennas  Prop. 
Vol.  59,  p.  1638,  2011).  Even-thought,  these  methods  have  offer  significant 
improvements,  there  are  still  cases  that  iterative  convergence  is  problematic 
or  simply  the  computational  overhead  is  too  large. 

The  talk  will  begin  by  exploring  the  landscape  of  the  state-of-the-art  DD 
methods  for  CEM.  In  the  main  part  of  the  talk  key  insights  about  the  root 
causes  of  the  robustness  and  scalability  problems  of  DD  methods  will  be  given, 
leading  the  way  to  the  development  of  a  new  class  of  preconditioners,  the 
locally  exact  algebraic  preconditioner  (LEAP)  and  the  Dual  Overlaping  Primal 
FETI  (FETI-DOP).  The  talk  will  close  with  results,  future  outlooks,  and 
extensions  to  outer  disciplines. 
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FEKO  (www.feko.info)  has  been  established  as  a  versatile  and  comprehensive  electromagnetic  simulation 
software  tool.  The  core  technology  is  the  method  of  moments  (MoM).  Over  the  past  two  decades  the 
MoM  has  been  hybridized  in  FEKO  with  other  methods,  including  the  Finite  Element  Method  (FEM), 
Physical  Optics  (PO),  Geometrical  Optics  (GO)  and  the  Uniform  Theory  of  Diffraction  (UTD). 
Hybridization  expands  the  range  of  application,  allowing  solution  of  problems  which  would  be  intractable 
with  any  individual  method.  Furthermore,  the  Multilevel  Fast  Multipole  Method  (MLFMM)  —  a  fast 
integral  equation  solution  technique  —  is  supported  in  the  MoM  and  hybrid  FEM/MoM  frameworks. 

This  presentation  will  give  an  overview  of  the  hybrid  FEM/MoM  capability  in  FEKO. 

The  hybrid  FEM/MoM  takes  the  full  coupling  between  metallic  structures  in  the  MoM  region  and 
dielectric  bodies  in  the  FEM  region  into  account.  The  strengths  of  the  two  techniques  can  thus  be 
exploited,  viz.  the  efficiency  of  the  MoM  for  modeling  open  boundary  radiating  structures  without 
discretization  of  the  environment,  and  the  efficiency  of  the  FEM  for  modeling  inhomogeneous  dielectric 
bodies.  The  MLFMM  is  an  acceleration  technique  for  the  MoM.  It  extends  the  capability  to  compute  full 
wave  solutions  for  electrically  large  structures,  and  in  particular  in  combination  with  the  FEM  for 
electrically  large  structures  together  with  complex  dielectrics  objects. 

Computational  aspects  regarding  the  FEM/MoM  and  FEM/MLFMM  include  the  following:  Iterative 
solvers  are  employed  together  with  sophisticated  preconditioners  for  an  efficient  computation.  A  good 
preconditioner  is  essential  to  ensure  convergence,  with  trade-offs  between  calculation  time,  memory 
usage  and  iterative  convergence  rate.  Parallel  processing  for  distributed  memory  environments  is 
supported.  Parallel  processing  achieves  faster  solution  times  and  by  employing  distributed  memory  larger 
problems  can  be  solved  on  less  expensive  computing  clusters. 

Electromagnetic  modeling  capabilities  for  the  FEM  include  the  following:  dielectric  and  magnetic 
materials,  perfect  conductor  surfaces,  complex  surface  properties,  e.g.  metallic  surfaces  exhibiting  the 
skin  effect,  modal  ports  of  arbitrary  cross  section,  impressed  filament  current  sources,  lumped  complex 
load  impedances.  The  FEM  and  MoM  can  be  decoupled  for  a  simplified  solution,  then  only  coupling 
from  MoM  to  FEM  is  taken  into  account  and  an  absorbing  boundary  condition  is  applied  to  the  FEM  and 
symmetry  may  be  employed  to  enhance  solver  efficiency. 

Typical  applications  include  antenna  simulations,  radiation  hazard  investigations  where  humans  interact 
with  RF  equipment,  the  modeling  of  waveguide  filters  with  dielectric  blocks  in  the  filter  and  the  modeling 
of  microstrip  patches  on  finite  substrates. 
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The  increasing  prevalence  of  handheld  mobile  devices  such  as  smartphones  and  tablet  computers  provides 
both  opportunities  and  challenges  to  antenna  design  engineers  nowadays,  because  the  designers  not  only 
need  to  have  antennas  cover  multi-band  and/or  multi-purpose  communications  as  standalone  elements, 
they  are  also  required  to  make  sure  the  antennas  perform  well  after  the  integration  into  the  devices.  Very 
often,  full  wave  analysis  needs  to  be  carried  out  to  study  the  mutual  couplings  between  multiple  antennas 
in  addition  to  their  radiation  characteristics,  considering  the  effects  not  only  from  handset  case  and  LCD, 
but  also  from  many  other  components,  such  as  battery,  speaker,  PCB,  etc. 

The  finite  element  method  (FEM)  is  intrinsically  suitable  for  such  problems  since  it  is  very  efficient  on 
modeling  3D  volumes  that  have  complex  geometrical  details  and  material  inhomogeneities.  However,  as 
FEM  formulation  doesn’t  incorporate  the  Sommerfeld  radiation  condition,  usually  the  radiation  condition 
in  FEM  implementations  is  applied  using  appropriate  absorbing  boundary  condition  (ABC),  or  perfectly 
matched  layer  (PML)  technique.  This  requires  the  domain  of  discretization  be  extended  far  from  the 
source  region.  The  integral  equation  method  such  as  Method  of  Moments  (MoM)  is  an  ideal  formulation 
for  open  boundary  radiating/scattering  structures.  As  a  result,  the  domain  of  discretization  can  be  kept  to  a 
minimum.  The  hybrid  FEM/MoM  exploits  the  benefits  from  both  techniques:  the  efficiency  of  the  MoM 
for  the  treatment  of  open  boundary  radiating  structures  because  no  discretization  of  three-dimensional  (3- 
D)  space  is  required,  and  the  efficiency  of  the  FEM  for  the  treatment  of  inhomogeneous  dielectric  bodies. 
The  approach  proves  to  be  especially  efficient  when  both  transmitter  and  receiver  need  to  be  included  in 
the  same  modeling  session.  It  is  also  a  preferred  method  for  bioelectromagnetic  problems  such  as 
radiation  hazard  with  specific  absorption  rate  (SAR)  calculations  (F.J.C.  Meyer,  etc.  “Human  exposure 
assessment  in  the  near  field  of  GSM  base-station  antennas  using  a  hybrid  finite  element/method  of 
moments  technique,”  IEEE  Trans.  Biomedical  Engineering,  vol.  50,  no.  2,  pp.  224-233,  Feb.  2003). 

In  this  talk,  various  examples  and  extensive  results  from  full  wave  simulations  with  hybrid  FEM/MoM 
for  handheld  mobile  devices  will  be  presented.  Performance  of  standalone  antennas  and  possible 
degradation  after  the  integration  will  be  examined.  Computational  resources  and  solution  accuracy  will  be 
discussed  as  well,  along  with  the  considerations  on  how  problem  complexity  and  material  inhomogeneity 
affect  the  computation  performance. 
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Abstract: 


The  finite  element  method  (FEM)  is  widely  used  in  computational  electromagnetics  since  it  is  very 
efficient  on  modeling  3D  volumes  that  have  complex  geometrical  details  and  material 
inhomogeneities.  For  radiation/scattering  problems,  usually  the  radiation  condition  in  FEM 
implementations  is  applied  using  appropriate  absorbing  boundary  condition  (ABC),  or  perfectly 
matched  layer  (PME)  technique.  This  requires  the  domain  of  discretization  be  extended  far  from 
the  source  region.  The  integral  equation  method  such  as  Method  of  Moments  (MoM)  is  an  ideal 
formulation  for  open  boundary  radiating/scattering  structures.  As  a  result,  the  domain  of 
discretization  can  be  kept  to  a  minimum.  Hybrid  FEM/MoM  technique  was  used  in  the  past  to 
predict  characteristics  of  cavity  backed  antennas  embedded  in  infinite  ground  plane.  In  most 
antenna-pattern  measurements,  the  infinite  ground  plane  is  simulated  by  a  relatively  large  finite 
ground  plane.  Although  this  simulation  does  not  appreciably  affect  the  input  admittance 
measurements,  diffractions  from  the  edges  of  the  ground  plane  modify  the  radiation  pattern 
considerably.  In  the  past,  the  geometrical  theory  of  diffraction  (GTD)  was  successfully  applied  to 
predict  the  changes  in  the  radiation  pattern  of  the  aperture  antennas  caused  by  a  finite  ground  plane. 
In  ref  1  and  2,  various  cavity  backed  apertures  (a)  open  coaxial  line1,  (b)  coaxial  cavity1  and  (c) 
cavity-backed  spiral  antenna2  were  analyzed  using  FEM/MoM/GTD  technique.  The  hybrid 
FEM/MoM  was  introduced  with  the  commercial  code  FEKO.  Simulations  are  carried  out  to 
validate  implementation  hybrid  FEM/MoM  technique  using  measured  data  in  ref  (2).  Some  of  the 
validation  results  were  presented  at  FEM  2006  workshop,  Though  we  have  seen  excellent 
agreement  between  measured  data  and  FEKO  simulations  for  cases  (a)  and  (b),  we  noticed 
discrepancy  between  FEKO  simulations  and  measured  data  for  case  (c).  In  this  paper  we  revisited 
this  problem  to  investigate  and  resolved  discrepancies  observed  earlier.  In  addition,  we  also  applied 
newly  introduced  FEM/MEFMM  technique  in  FEKO  to  the  above  problems  for  validation  using 
the  experimental  data  in  ref  (1  &  2).  Numerical  results  will  be  presented  at  the  FEM  2012 
workshop. 


1  C.J.Reddy,  M.D.Deshpande,  C.R.Cockrell  and  F.B.Beck,  “Radiation  characteristics  of  cavity  backed 
aperture  antennas  in  finite  ground  plane  using  hybrid  FEM/MoM  technique  and  Geometrical  Theory  of 
Diffraction,”  IEEE  Transactions  on  Antennas  and  Propagation,  Vol.44,  No. 10,  pp. 1327-1333,  October  1996. 

2  C.J.Reddy,  M.D.Deshpande,  D.T.Fralick,  C.R.Cockrell  and  F.B.Beck,  “A  combined  FEM/MoM/GTD 
technique  to  analyze  elliptically  polarized  cavity-backed  antennas  with  finite  ground  plane,”  NASA 
Technical  Paper  3618,  November  1996. 
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Design  of  high  power  filters  and  multiplexers  for  modern  communication  systems  involves  full  wave 
optimization  based  on  extensive  electromagnetic  simulations  of  complex  3D  structures.  These 
simulations  are  often  carried  out  using  finite  element  method  (FEM)  due  its  ability  to  handle 
arbitrary  shapes.  However,  the  numerical  cost  associated  with  FEM  analysis  is  high  as  one  has  to 
evaluate  filter’s  response  at  several  frequencies  and  optimize  the  design  to  fulfil  the  specifications, 
while  modifying  filter  geometry. 

A  classic  approach  for  computation  of  scattering/impedance  parameters  for  a  structure  using  FEM 
leads  to  a  linear  system  A-x  —  b  that  has  to  be  solved  anew  at  each  frequency  point  for  a  collection 
of  excitations  at  the  feeding  ports  (for  the  TV-port  device  one  gets  N  right  hand  sides  at  each 
frequency).  Since  matrix  A  is  usually  large,  sparse  and  frequency-dependent,  the  numeric  cost  of 
solution  procedure  (using  either  direct  or  iterative  methods)  is  high,  especially  if  response  at  tens 
or  hundreds  frequency  points  needs  to  be  computed. 

An  alterative  approach  is  to  replace  the  repeated  solutions  of  linear  problems  with  a  defined  set  of 
eigenproblems.  The  new  formulation  binds  the  eigenvalues  of  the  predefined  eigenproblems  with 
zeros  and  poles  of  the  rational  representation  of  impedance  parameters  ^j(s)  =  of  the  ana¬ 
lyzed  structure.  Additionally,  it  is  shown  that  the  resulting  eigenproblems  can  be  efficiently  solved 
using  iterative  eigensolvers,  such  as  shift-and-invert  Lanczos  methods.  The  proposed  formulation 
is  particularly  well  suited  for  analysis  of  high  frequency  filtering  devices,  since  their  responses  are 
commonly  described  by  zeros  and  poles  of  the  transfer  function.  This  is  because  the  extracted  zeros 
and  poles  can  be  directly  used  in  optimization  of  microwave  filters  using  a  special  cost  function 
involving  sets  of  reference  zeros  and  poles  of  ideal  filter  prototype  (  P.  Kozakowski,  M.  Mrozowski, 
’’Automated  CAD  of  coupled  resonator  filters,”  IEEE  Microwave  and  Wireless  Components  Let¬ 
ters,  vol.12,  no.  12,  pp.  470-  472,  Dec  2002).  The  applications  of  the  new  methodology  based  on 
eigenvalue  computations  for  optimization  of  microwave  filters  using  FEM  will  be  presented  at  the 
workshop. 
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Fast  frequency  sweep  (FFS)  is  a  common  term  for  highly  efficient  solution  strategies  for  the  broad¬ 
band  finite- element  (FE)  analysis  of  passive  microwave  structures.  However,  many  applications  do 
not  only  require  the  frequency  response  of  certain  system  outputs  but  also  their  sensitivities  w.r.t. 
variations  in  parameters  such  as  geometry  or  material  properties.  In  case  of  the  original  FE  system, 
sensitivities  at  a  given  solution  frequency  are  computed  very  efficiently  by  means  of  adjoint  variables 
techniques:  a  number  of  forward-back  substitutions  suffice.  It  is  therefore  very  appealing  to  apply 
the  same  method  to  the  reduced-order  systems  of  FFS’s  [1],  [2].  Even  though  this  approach  has 
been  used  very  successfully  in  the  context  of  FE’s  [3],  its  convergence  is  far  from  obvious,  for  the 
following  reasons:  First,  sensitivities  are  no  longer  given  by  first-order  derivatives  about  expansions 
points  at  which  the  exact  solution  of  the  FE  system  is  known.  Second,  the  reduced-order  models 
of  FFS’s  are  constructed  with  only  the  frequency  in  mind.  Since  the  trial  and  testing  spaces  are  of 
very  low  order,  they  may  hold  little  information  about  derivatives  w.r.t.  to  some  other  parameter. 

The  first  part  of  the  talk  gives  a  formal  convergence  analysis  in  the  general  context  of  projection- 
based  MOR,  which  covers  both  single-  and  multi-point  methods.  It  will  be  proven  theoretically 
and  demonstrated  numerically  that,  under  the  precondition  that  two-sided  projections  are  used, 
i.e.  that  both  excitation  vectors  and  output  functionals  are  taken  into  account,  sensitivities  of 
these  outputs  do  converge,  at  a  lower  rate  though.  However,  sensitivities  of  the  (generalized)  state 
variables  do  not  converge.  The  latter  deficiency  is  cumbersome  in  those  cases  when  outputs  are 
constructed  in  a  post-processing  step.  Typical  applications  include  the  a  posteriori  definition  of 
field  probes,  Huygens  surfaces,  or  similar.  The  second  part  of  the  talk  introduces  an  improved 
FFS  approach  that  guarantees  convergence  of  the  sensitivities  of  the  (generalized)  system  state  and 
hence  of  any  output  defined  in  a  post-processing  step.  We  will  present  a  general  framework  for  the 
construction  of  trial  and  testing  spaces  which  also  covers  recent  works  in  circuit  theory  [4],  and 
give  numerical  examples  that  illustrate  our  findings. 
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This  work  introduces  a  conformal  finite  difference  method  (CFDM)  to  accurately  and  efficiently  study 
resonances  in  electromagnetic  or  radio  frequency  (RF)  structures.  For  benchmark,  we  calculate  the 
resonances  of  a  circular  cavity  using  the  CFDM  and  compare  the  results  with  those  from  the  finite 
element  method  (FEM).  The  accuracy  of  the  CFDM  is  measured  by  comparing  with  calculations  based  on 
the  finite  element  method.  The  results  show  that  an  accuracy  of  99.4%  can  be  achieved  by  using  only 
10,000  mesh  points  with  the  Dey-Mittra  algorithm  as  implemented  in  the  CFDM.  By  comparison,  a  mesh 
number  of  250,000  is  needed  to  preserve  99%  accuracy  using  a  staircased  FDM.  This  suggests  one  can 
more  efficiently  and  accurately  study  RF  structures  using  the  CFDM  than  a  conventional  FDM.  For 
illustration,  the  dispersion  relation  of  an  A6  relativistic  magnetron  has  also  been  determined  and  a 
comparison  with  finite  element  method  has  been  done. 
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Motivation 

With  the  continuous  advancements  in  semiconductor  technology,  IC  circuit’s  operating  speed,  routing  density, 
and  routing  complexity  have  increased  rapidly  while  the  operating  voltage  and  voltage  swing  have  decreased 
sharply.  Unwanted  electromagnetic  behaviors,  such  as  signal  reflection,  cross-talk,  signal  noise,  and 
electromagnetic  interference  impact  significantly  the  performance  of  the  circuitry  functionalities  and  lead  to  the 
misjudgment  of  the  digital  signal’s  level  and  timing.  The  desired  and  unwanted  electromagnetic  behavior  over 
interconnects  among  different  subsystems,  functional  modules  on  a  PCB  or  in  a  packaging  must  be  analyzed  and 
characterized  accurately  and  efficiently  to  ensure  a  reliable  system  design. 

The  ever-increasing  speed  and  the  interconnect  complexity  lead  to  a  wavelength  comparable  to  the 
interconnect  size  which  necessitates  full- wave  electromagnetic  analysis.  However,  such  a  full- wave  analysis  is  in 
general  highly  computationally  demanding  in  terms  of  both  memory  and  computation  time.  Most  of  the  full-wave 
analysis  algorithms  typically  require  large-scale  matrix  equation  solving  which  consumes  the  most  memory  and 
CPU  usage.  It  has  been  one  of  the  most  crucial  challenges  to  effectively  solve  the  matrix  equations  with  hundreds 
of  thousands  and  even  millions  of  unknowns  to  make  the  full- wave  characterization  and  simulation  practical.  One 
way  to  reduce  the  complexity  is  domain  decomposition.  In  this  presentation,  we  discuss  two  decomposition 
mechanism  implemented  Cadence  Allegro  product  line  with  finite  element  method.  One  is  layout-aware 
decomposition  and  the  other  is  EM  solver  specific  decomposition. 

Layout  Level  Domain  decomposition 

Typically  in  a  PCB  design,  an  initial  layout  exists  which  in  parallel  requires  interactively  SI  analysis.  Figure  1 
shows  a  layout  where  a  signal  channel  is  selected  for  characterization  and  SI  analysis.  Look  into  the  highlighted 
signal  channel,  clearly  it  would  be  resource  misuse  to  solve  the  entire  channel  as  one  3D  structure  when  it  can  be 
decomposed  into  sub-domains  while  on  the  decomposition  interfaces  a  few  propagation  modes  can  be  well 
defined  and  be  easily  found.  This  layout  level  of  decomposition  can  be  applied  to  complicated  SiP  or  SOP 
designs  in  which  systems  are  contained  in  the  other  system  and  are  implemented  in  one  or  multiple  stacked 
packages. 


Figure  1  Layout  in  Cadence  Allegro  editor  Figure2  Sub-domain  for  EM  Solver  to  apply 

Domain  decomposition  Within  EM  Solver 

For  each  decomposed  sub-domains,  computational  electromagnetic  methods  will  be  applied  to  extract 
necessary  electrical  parameters,  either  in  FD  or  in  TD.  Since  each  so  divided  sub-domain  can  still  be  large  and 
complicated,  respected  to  today’s  computation  software  and  hardware,  further  domain  decomposition  needs  to  be 
performed  in  the  EM  computation  level.  For  this  purpose,  the  finite  element  method  is  highly  desired  and 
efficient.  Cadence  and  UIUC  are  jointly  working  actively  on  these  levels  of  decomposition.  Those  sub-domains 
characterization  or  simulation  is  performed  parallel  on  multi-core  and  even  on  distributed  computation  systems. 

SUB-DOMAIN  INTEGRATION 

The  decomposition  is  performed  in  such  a  way  that  each  sub-domain  is  constructed  in  SPICE  sub-circuit 
format  and  all  the  I/Os  of  those  sub-circuit  are  clearly  defined  in  the  layout  level  structure  extraction.  The  system 
level  simulation  can  be  performed  naturally.  Cadence  and  UIUC  has  developed  full-wave  spice  to  address  the 
complicated  EM  behavior  on  SiP  design  which  includes  EMI  and  EMC  prediction. 
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Circuit  simulation  is  an  increasingly  indispensable  tool  for  the  design  of  integrated  circuits  (ICs).  The 
most  prominent  circuit  simulation  programs  are  SPICE  (Simulation  Program  with  Integrated  Circuit 
Emphasis)  and  its  derivatives.  SPICE  is  highly  capable  of  simulating  active  devices.  In  the  early  years 
when  interconnects  could  be  modeled  simply  as  lumped  capacitors,  the  linear  network  is  a  very  small 
component  of  an  IC.  The  exponentially  increased  complexity  of  ICs,  however,  has  made  circuit 
simulation  increasingly  challenging. 

In  the  prevailing  circuit  simulation  paradigm,  even  though  the  simulation  of  a  circuit  is  the  end  goal,  the 
state-of-the-art  circuit  simulation  flow  involves  first  the  construction  of  an  RLC-based,  Transmission- 
line -based,  or  S-parameters-based  representation  of  the  linear  network;  this  step  is  called  extraction.  The 
extracted  model  of  the  linear  network  is  then  stitched  with  the  nonlinear  devices  for  the  simulation  of  the 
entire  circuit.  The  sheer  dominance  of  the  linear  network  over  the  nonlinear  network  in  today’s  and  future 
generation  integrated  circuits  poses  challenges  to  the  prevailing  simulation  paradigm  in  both  extraction 
and  simulation. 

In  this  work,  we  develop  an  extraction- free  circuit  simulator  that  overcomes  the  challenges  facing  existing 
SPICE-based  circuit  simulators.  It  bypasses  the  step  of  extraction,  while  retaining  an  RLC-based 
perspective  that  allows  a  designer  to  continue  to  exercise  his  or  her  circuit  intuition.  The  underlying 
model  for  simulation  is  inherently  passive  and  stable.  Moreover,  it  captures  the  interaction  between  the 
nonlinear  circuit  and  the  linear  network,  with  the  nonlinear-linear  coupling  rigorously  taken  into 
consideration.  Perhaps  more  important,  it  allows  for  the  simulation  of  a  circuit  of  0(N)  size  including 
both  nonlinear  devices  and  the  linear  network  in  0(N)  complexity,  i.e.,  optimal  complexity.  As  an  added 
bonus,  the  proposed  circuit  simulator  achieves  electromagnetic-physics  based  accuracy.  Application  to 
very  large-scale  on-chip  circuits  involving  close  to  1  million  CMOS  transistors  and  interconnects  having 
hundreds  of  millions  of  unknowns  has  demonstrated  the  superior  performance  of  the  proposed  extraction- 
free  circuit  simulator. 
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We  consider  a  parameterized  multiple  electromagnetic  scattering  wave  propagation  model  in  three 
dimensions.  The  parameters  in  the  model  describe  the  location,  orientation,  size,  shape,  and  number 
of  scattering  particles  as  well  as  properties  of  the  input  source  field  such  as  the  frequency,  polariza¬ 
tion,  and  incident  direction.  The  need  for  fast  and  efficient  (online)  simulation  of  the  interacting 
scattered  fields  under  parametric  variation  of  the  multiple  particle  surface  scattering  configuration 
is  fundamental  to  several  applications  for  design,  detection,  or  uncertainty  quantification. 


For  such  dynamic  parameterized  multiple  scattering  models,  the  standard  discretization  procedures 
are  prohibitively  expensive  due  to  the  computational  cost  associated  with  solving  the  full  model 
for  each  online  parameter  choice.  In  this  work,  we  propose  an  iterative  offline/online  reduced  basis 
computer  model  reduction  approach  for  a  boundary  element  method  to  simulate  a  parameterized 
system  of  surface  integral  equations  reformulation  of  the  multiple  particle  wave  propagation  model. 


The  approach  includes  (i)  a  greedy  algorithm  based  computationally  intensive  offline  procedure  to 
create  a  selection  of  a  set  of  a  snapshot  parameters  and  the  construction  of  an  associated  reduced 
boundary  element  basis  for  each  reference  scatterer  and  (ii)  an  inexpensive  online  algorithm  to 
generate  the  surface  current  and  scattered  field  of  the  parameterized  multiple  wave  propagation 
model  for  any  choice  of  parameters  within  the  parameter  domains  used  in  the  offline  procedure. 
Comparison  of  our  numerical  results  with  experimentally  measured  results  for  some  benchmark 
configuration  demonstrate  the  power  of  our  method  to  rapidly  simulate  the  interaction  of  scattered 
wave  fields  under  parametric  variation  of  the  overall  multiple  particle  configuration. 
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